Spacecraft Radiation Mapping

[CharlesPlusC]

Hello Mitsuba Community,

I am a PhD student utilising Mitsuba for my research on solar radiation pressure (SRP) modelling, specifically for spacecraft. My academic foundation is in astrodynamics, not rendering; hence, I appreciate your patience with any naive questions or oversights.

Objective:
I aim to generate a comprehensive map depicting both the incoming and outgoing direction and magnitude of solar radiation across a spacecraft's surface. This will enable me to determine the cumulative force exerted by solar radiation.

Current Methodology:
1.Scene Definition: Construct a Mitsuba scene incorporating:

• A directional light source (type: spectrum, value: 1367w/m^2).
• An .obj file showcasing the target object, currently utilising a preset material property.
• The ptracer integrator.

2.Surface Interaction Points: With a camera, I accumulate multiple surface interaction points on the object, saving their shader frames and positions. I rotate the camera around the object, continually documenting these points. Subsequently, I resample these points for even distribution (I specify the spacing myself) on the object's surface using KDTrees and save the s, t, n vectors for each.

For ray generation:

ray_origin_local = mi.Vector3f(x, y, 0)
ray_origin = mi.Frame3f(cam_dir).to_world(ray_origin_local) + cam_origin
camera_rays = mi.Ray3f(o=ray_origin, d=cam_dir)
si = scene.ray_intersect(camera_rays)

3.Hemisphere of Sensors: I then generate a unit hemisphere of orthographic sensors (where the number of sensors in this sensor hemisphere is defined manually).

4.Sensor Application: I apply this hemisphere of sensors to each of the surface interaction points in the list (from step 3). So if my object ends up having 700 surface interaction points, and I decide to make 10 hemisphere partitions per point, I will have 7000 sensors. Since each shader frame has a local s,t,n reference frame, I rotate the hemisphere to be in the correct orientation for each of these points.

5.Rendering: Then I render each of these sensors, and take the average value of all the pixels in the rendered image (I wonder if this is not part of my problem).

6.Vector Averaging: Now for each point on my object I have radiation 10 vectors, each pointing in the direction that the sensor was pointing and scaled by the average strength of the pixels in the image rendered. I average the 10 vectors at each location so that each surface interaction point on the surface of my object now has an associated magnitude and strength of radiation.

Issue:
While the relative magnitude of the "radiation maps" align with my expectations for several test cases, the absolute magnitude is inaccurate and appears to vary with the number of interaction points and the number of hemisphere cameras.

I believe that this may be in part because the sensors at each location on the surface of the object are not rendering the light directly coming from the emitter, but rather only the light that is reflected from the surface of the object. I have come to this conclusion after inspecting the individual renders that are pointing to the light source: they are blank (black).

I've perused the documentation, but I'm still uncertain about how to configure the sensors to capture this direct light or if there is a better of doing this altogether... Despite experimenting with various integrators and light sources, I've had no success. I'm contemplating if a distinct sensor type might be the solution. The irradiance meter seemed to be a good match for my needs but I was not able to make it work on any .obj.

Please find attached sample outputs of my method on a unit sphere:

I trust this summary clearly conveys my efforts. While I recognise that a thorough review of my approach might be demanding, I'd deeply appreciate any insights or suggestions you can offer.

Thank you for your time, and for offering this exceptional tool to the community.

Best regards,
Charles

[njroussel]

Hi @CharlesPlusC

This is definitely not something I was expecting 😄 -- it's cool to see Mitsuba being used in such a context.

I do need a few clarifications:

1. First of all, you mention that you are using the ptracer integrator. Is this integrator used in step 5? If so, I'm struggling to see how this could possibly work. The orthographic plugin is missing a few methods to "connect" light paths to it. It should crash with an error message about some method not being implemented. Unless you've made some changes of your own to either plugin.
2. This is a follow-up question. Are you using a stock build of Mitsuba? Especially, I'd need to know if you made any changes to the sensors, integrators, or emitters you are using.

not rendering the light directly coming from the emitter

This is somewhat expected. Your scene is built with a few hard constraints: the orthographic sensor, and directional emitter can respectively only receive and emit light from a single direction. Without a surface to bounce-off of, these directions will never match and hence a light path cannot be constructed from the emitter to the sensor or the other way around.
Of course, there is an edge case where both are perfectly aligned. I would need to double-check the code, but I am practically sure this case is not handled correctly as we assume that it is statistically impossible to find such a light path.

The irradiance meter seemed to be a good match for my needs but I was not able to make it work on any .obj.

The irradiancemeter measures the incident radiance across the entire shape it is attached to. You would need to define discrete areas on your shape/spacecraft, or else you would be capturing the sum over the entire shape.
In addition, you would not be able to capture the directionality.

[CharlesPlusC]

Hi @njroussel,

Thank you so much for your swift response! 😃

you mention that you are using the ptracer integrator

My apologies for the oversight. While I did employ ptracer in other steps, I actually used the path integrator for Step 5. I understand now that I had mixed up "path tracer" and "particle tracer."

Are you using a stock build of Mitsuba?

Indeed, I'm working with the stock version of Mitsuba. The only alteration I've made is setting the variant to llvm_ad_rgb.

Given the results I've observed, would it be accurate to describe my current output as a map showcasing reflected radiation, without accounting for the direct incoming radiation?

You would need to define discrete areas on your shape/spacecraft

Interesting suggestion! Thank you. For capturing the directionality of incoming radiation, do you have any recommended approaches? I've considered "wrapping" the spacecraft in a transparent material. Perhaps by creating a duplicate of its geometry but adjusting the material properties to be "inert" during light interactions. My hypothesis is that if light travels through this transparent layer (and if this qualifies as a "bounce"), it might then be more likely to intersect with a sensor.

Thank you once more for your insights and assistance.

Warm regards,
Charles

[njroussel]

Hi again,

Sorry it took some time to come back to this thread.

Given the results I've observed, would it be accurate to describe my current output as a map showcasing reflected radiation, without accounting for the direct incoming radiation?

Yes, but not quite. One last problem is that you're not measuring the reflected radiation of a single point. The hemisphere of sensors can obviously see a bit around the actual point of interest and you will be capturing that radiance too. (You can minimize this problem by having the hemisphere extremely close to the point).

If you want numerical correctness, my recommendation would be to implement a hemispherical sensor as a custom plugin (tutorial). Actually, a similar ideas was discussed in this thread with an implementation. In that use case, a spherical sensor with the origin of the rays at the center of the sphere is needed. In your case, you'd want a hemisphere or ray origins, all pointing to the same direction. This should be fairly easy to implement.
One last point is that with a hemispherical sensor, your output is still a 2-dimensional rectangular film/tensor when rendering. You need to convert this 2D planar information to your hemisphere. The parameterization you defined in the custom sensor to convert a random 2D sample to a direction on the hemisphere will be exactly this "conversion". This will allow you to compute the average reflected direction.

[CharlesPlusC]

Hi @njroussel,

Thank you for your reply. There’s absolutely no need to apologize for the delay; I'm just very grateful to have the benefit of your expertise!

my recommendation would be to implement a hemispherical sensor as a custom plugin

OK, I will give this a shot in the coming weeks. I appreciate the suggestion and the link to the relevant thread.

Once again, a big thank you for your input. I will be sure to credit you in any work that comes out of this.

[CharlesPlusC]

Hi @njroussel, 👋 😄

I hope you are well!

I'm back with a few more queries, and I must apologize for the length of this post. Your insights thus far have been immensely valuable, and I completely understand if your time constraints don't permit a detailed response to these follow-up questions.
Regardless, I thought it'd be beneficial to outline my thought process for the community's reference, hoping it might spark further discussion or guidance from anyone willing to contribute.

For reference, here is the current state of my hemispherical sensor plugin:

def square_to_uniform_hemisphere(sample, scale_factor=1.0):
    phi = 2 * sample.x * dr.pi

    # Use the second sample dimension to get sin(theta) and cos(theta)
    cos_theta = sample.y
    sin_theta = dr.sqrt(1.0 - cos_theta * cos_theta)

    # Convert spherical to Cartesian coordinates
    x = dr.cos(phi) * sin_theta 
    y = cos_theta
    z = -dr.sin(phi) * sin_theta

    return scale_factor * mi.Vector3f(x, y, z)
    
class HemisphericalCamera(mi.Sensor): 
    """Defines a hemispherical sensor with inward-pointing rays."""
    ray_origins = []
    ray_directions = []
    RTN = np.empty((0, 3, 3), float)

    def __init__(self, props=mi.Properties()):
        super().__init__(props)
        #scale the size of the sensor relative to world units
        self.scale_factor = props.get('scale_factor', 1) 
        # Add the hemisphere_center to the camera properties 
        self.hemisphere_center = props.get('hemisphere_center', mi.Vector3f(0, 0, 0))
        # Compute the RTN vectors for this sensor instance (for visualization)
        self.compute_and_store_RTN()

    def compute_and_store_RTN(self):
        # Assuming the default 'up' direction and 'north' for the sensor when it's created
        up = np.array([0.0, 1.0, 0.0])  # Y-up convention
        north = np.array([0.0, 0.0, 1.0])  # Adjusted "north" towards positive Z since we're using a left-handed system in Mitsuba

        # Calculate the RTN vectors based on 'up' and 'north'
        r = up  # Radial vector pointing outwards from the dome's base
        t = np.cross(north, up)  # Changed the order for left-handed system
        t /= np.linalg.norm(t)  # Normalize the 't' vector
        n = np.cross(r, t)  # Normal, points towards the dome center, perpendicular to 'r' and 't'

        # Convert the RTN vectors from NumPy arrays back to Mitsuba's Vector3f objects
        r_mitsuba = mi.Vector3f(r[0], r[1], r[2])
        t_mitsuba = mi.Vector3f(t[0], t[1], t[2])
        n_mitsuba = mi.Vector3f(n[0], n[1], n[2])

        #convert to world coordinates
        r_vec = self.world_transform() @ r_mitsuba
        t_vec = self.world_transform() @ t_mitsuba
        n_vec = self.world_transform() @ n_mitsuba

        rtn_array = np.array([r_vec, t_vec, n_vec], dtype=float).reshape(1,3,3)

        # Append the RTN vectors for this instance to the class-level list
        HemisphericalCamera.RTN = np.concatenate((HemisphericalCamera.RTN, rtn_array), axis=0)

    def sample_ray(self, time, wavelength_sample, position_sample, aperture_sample, active=True):
        #this function is called during the path tracing integrator call
        
        wavelengths, wav_weight = self.sample_wavelengths(dr.zeros(mi.SurfaceInteraction3f),
                                                          wavelength_sample, active)

        center = self.world_transform().translation()
        local_o = square_to_uniform_hemisphere(position_sample, self.scale_factor)
        o = self.hemisphere_center  # Using the hemisphere_center property to adjust the ray's origin

        o += self.world_transform() @ local_o
        d = center - o
        
        # Append the ray's origin and direction to the class-level lists for visualization
        self.ray_origins.append(o)
        self.ray_directions.append(d)

        return mi.Ray3f(o, d, time, wavelengths), wav_weight

    def sample_ray_differential(self, time, wavelength_sample, position_sample, aperture_sample, active=True):
        ray, weight = self.sample_ray(time, wavelength_sample, position_sample, aperture_sample, active)
        return mi.RayDifferential3f(ray), weight

    def sample_direction(self, it, sample, active=True):
        # this function is called during the particle tracing integrator call
        # this is not useful though because the rays are not sampled according to the "square_to_uniform_hemisphere" function
        trafo = self.world_transform()
        ref_p = trafo.inverse() @ it.p
        d = mi.Vector3f(ref_p)
        dist = dr.norm(d)
        inv_dist = 1.0 / dist
        d *= inv_dist
        resolution = self.film().crop_size()

        ds = dr.zeros(mi.DirectionSample3f)

        ds.uv = mi.Point2f(dr.atan2(d.x, -d.z) * dr.inv_two_pi, dr.safe_acos(d.y) * dr.inv_pi / 2.0)
        ds.uv.x -= dr.floor(ds.uv.x)
        ds.uv *= resolution
        sin_theta = dr.safe_sqrt(1 - d.y * d.y)

        ds.p = trafo.translation()
        ds.d = (ds.p - it.p) * inv_dist
        ds.dist = dist
        ds.pdf = dr.select(active, 1.0, 0.0)

        weight = (1 / (2 * dr.pi * dr.pi * dr.maximum(sin_theta, dr.epsilon(mi.Float)))) * dr.sqr(inv_dist)
        
        return ds, mi.Spectrum(weight)

mi.register_sensor("hemispherical", lambda props: HemisphericalCamera(props))

Good news:

• The renders make sense when inspected individually (i.e. shadows and light where I expect them, magnitude of values increasing on surfaces receiving more light).

Bad news:

• The magnitude of the values from the sensor pixels vary drastically depending on the integrator, BSDF, and ray "lengths" (more on this below) and I do not fully understand why.

BSDF/Integrator Issues

Having reviewed some relevant posts, notably (#376) and (#365), it seemed ptracer should be the correct choice for what I am trying to do.

My initial attempt at deciding what integrator to use involved setting up a render with a path integrator, a directional emitter to mimic solar radiation in space:

    'emitter': {
        'type': 'directional',
        'direction': [1.0, 0.0, 0.0],
        'irradiance': {
            'type': 'rgb',
            'value': 1367,
    },
    },

hdrfilm film, independent sampler, and rendering a diffuse sphere with a perspective sensor:

    'sphere' : {
        'type': 'sphere',
        'radius': 1.0,
        'bsdf': { 
            'type': 'diffuse',
            'reflectance': {
                'type': 'rgb',
                'value': [1, 1, 1]
            }
        },
    },
    'sensor': {
        'type': 'perspective',
        'to_world': mi.ScalarTransform4f.look_at(origin=[-5, 0, 0],
                                                 target=[0, 0, 0],
                                                 up=[0, 0, 1]),

Querying the pixel at the center of the image I get a value of 435.12. This is very close to 1/3.14*(1367) which matches up with the analytical solution one would expect for the normal component of radiation reflected by a diffuse material.

Switching to a highly specular BSDF, however, results in physically implausible values.

    "RoughBSDF": {
        "type": "roughconductor",
        "material": "none",
        "alpha": 0.08,
        "distribution": "ggx",
    },
    'sphere': {
        'type': 'sphere',
        'radius': 1,
        "bsdf": {
            "type": "ref",
            "id": "RoughBSDF",
                },
}

With alpha=0.08, the readings skyrocket to 15088.054 and 13111.886 for the path and ptracer integrators, respectively, despite an irradiance = 1367 W/m2. Adjusting alpha = 0.28 results in a more plausible ~1374.165, but an alpha = 0.99 plummets the value to 99. This is unexpected as I would have thought roughconductor with a high alpha would mimic the diffuse scenario more closely.

I'm at a loss interpreting these results. Are any of these readings physically accurate? Does the 'path' integrator have an "upper" limit of 0.28 and lower limit of 0.9 for rendering specular and diffuse surfaces? Is such a calibration valid at all?

In addition to returning very similar values to the path integrator above (thus not solving the problems associated with directional light on a highly specular surface), my understanding is that ptracer is also incompatible with my hemispherical sensor since ptracer doesn't sample along predefined rays (such as defined in my hemisphere plugin) but rather starts from the light source attempting to connect with the object/sensor. I did try ptracer with my hemisphere plugin anyways but it lead to weird looking renders, hence I am sticking to path for now.

Hemisphere Issues

This is what one instance of my current sensor looks like (only 1/10000 rays plotted here) when positioned on a sphere directly opposite the incoming directional light source:

Placing 200 sensors around a unit sphere with seems to model the physical behaviour well enough within a certain range of alpha values [0.28-0.9], but anything outside these bounds starts to lose physical sense...

Low alpha (0.28) Highest mean sensor pixel value ~1400

High Alpha(0.85) Highest mean sensor pixel value ~420

ts/66725307/49794064-bccb-4fac-b4be-3f266684416a)

In my sensor plugin I have defined 2 "extra" parameters:

• scale_factor scales the hemisphere's radius, essentially dictating the rays' starting distance. This parameter was introduced to prevent rays from starting within another nearby shape in scenarios with more cluttered geometry.

• hemisphere_center serves as an offset, positioning the hemisphere above the sampled surface. This was added as I anticipated complications if the hemisphere's center was at the same position as the geometry, though this concern seems unfounded.

To my confusion, altering scale_factor modifies the recorded values. For instance, a scale_factor of 0.2 results in an average pixel value of 4837, while a scale_factor of 1 reduces it to 276. A scale_factor of 0.5 yields values closest to the anticipated (~1380 for the sensor directly opposite to the the incident radiation) so I have stuck with this for now but this is purely empirical and uninformed.
I'm puzzled by how the ray lengths could influence the recorded values, as my understanding was that the rays merely needed a non-zero length, since what is fundamentally captured is the light post-ray-surface interaction...

As always I'm genuinely grateful for any insight you can offer on this matter. I also want to express my commitment to appropriately and honestly acknowledging your contributions in this work. Should my work here lead to any future publications or presentations, I will ensure credit is given in a manner that accurately reflects your intellectual input and guidance, as per your wishes and academic norms.

Thank you once again for your time and expertise!

Charles

[njroussel]

Hi @CharlesPlusC 😄

---

With alpha=0.08, the readings skyrocket to 15088.054 and 13111.886 for the path and ptracer integrators, respectively, despite an irradiance = 1367 W/m2. Adjusting alpha = 0.28 results in a more plausible ~1374.165, but an alpha = 0.99 plummets the value to 99. This is unexpected as I would have thought roughconductor with a high alpha would mimic the diffuse scenario more closely.

I'm at a loss interpreting these results. Are any of these readings physically accurate? Does the 'path' integrator have an "upper" limit of 0.28 and lower limit of 0.9 for rendering specular and diffuse surfaces? Is such a calibration valid at all?

Low roughness / highly specular materials will be very noisy, you'll need a pretty high sample count to get a good estimate of the radiance.
As for high roughness materials, the model implement in Mitsuba is not perfectly energy conserving. There exist a line of work (this paper for example) which properly simulates multiple scattering events on the micro-scale geometry and will not lose energy. So, you are somewhat right - there is a range of alpha values that behaves "correctly". The rougher the material, the higher the energy loss.
Note: the diffuse BSDF is energy conserving, you can use it to validate your setup as you're doing currently.

---

In addition to returning very similar values to the path integrator above (thus not solving the problems associated with directional light on a highly specular surface), my understanding is that ptracer is also incompatible with my hemispherical sensor since ptracer doesn't sample along predefined rays (such as defined in my hemisphere plugin) but rather starts from the light source attempting to connect with the object/sensor. I did try ptracer with my hemisphere plugin anyways but it lead to weird looking renders, hence I am sticking to path for now.

Yes, I don't think there is a need for anything else than the path integrator for now.

---

To my confusion, altering scale_factor modifies the recorded values. For instance, a scale_factor of 0.2 results in an average pixel value of 4837, while a scale_factor of 1 reduces it to 276. A scale_factor of 0.5 yields values closest to the anticipated (~1380 for the sensor directly opposite to the the incident radiation) so I have stuck with this for now but this is purely empirical and uninformed.
I'm puzzled by how the ray lengths could influence the recorded values, as my understanding was that the rays merely needed a non-zero length, since what is fundamentally captured is the light post-ray-surface interaction...

I'm puzzled by this. Admittedly, I haven't looked at your sensor implementation in detail (by the way, thank you for sharing it publicly! 🎉). My immediate guess is that there is some term missing to account for the pixel's surface area. Or that maybe when the ray origins are too far from the surface, the intersection point slightly moves because of floating point errors (this seems less likely).
You can check that the actual ray contributions (Integrator::sample()) remain identical for two different scale_factor values, when using the same set of random numbers across both runs. You'll want to do this with a fairly low sample count and it will require adding print statement in the source code. If both values are identical, then it should just be a matter of following the source code from that point on to see when they diverge in values. If they aren't the same, and the rays are indeed the same, then something weird is happening inside the integrator's code.

[CharlesPlusC]

Hello @njroussel, 👋 🤠

As per your recommendation I developed a scaling function for the pixels and am pleased to share that the values I am now obtaining are making perfect physical sense across various geometries when using the diffuse BSDF 🎉 🥳

I have one more question that I hope you might be able to shed some light on. Could you please tell me if there is a list of the BSDFs in Mitsuba that do adhere to the principles of energy conservation? I have looked in the docs but was not able to find a clear list...

Looking ahead, provided I can find suitable BSDFs that approximate spacecraft material properties, my intention is to benchmark this method against state-of-the-art spacecraft radiation pressure mapping methods. All code generated in this process will be made available openly. I am keen on ensuring that your contributions to this endeavour are properly acknowledged in a manner that suits you best so I plan to reconnect once I have more results in hand to discuss your preferred mode of acknowledgment.

Once again, thank you for sharing your time and expertise!

Best,
Charles

[njroussel]

As per your recommendation I developed a scaling function for the pixels and am pleased to share that the values I am now obtaining are making perfect physical sense across various geometries when using the diffuse BSDF 🎉 🥳

Great 🎉 !

Could you please tell me if there is a list of the BSDFs in Mitsuba that do adhere to the principles of energy conservation?

We unfortunately don't have this. In short, any material that relies on microfacet theory (matierals which have a roughness parameter typically) are not energy-conserving. I think that leaves only the diffuse, dielectric and conductor.
Note that the exact energy loss is usually known and can be compensated for. Sometimes, this loss of energy is re-injected as a diffuse component to the BSDF. (See section 5.3 and figure 4 of this paper for example: http://www.cs.cornell.edu/projects/layered-sg14/layered.pdf)

[CharlesPlusC]

Perfect. Thanks very much again for your help with all this!

Best,
Charles