Automatically Generating struct Descriptions

It is often convenient and efficient to write raw data to files. Converting numbers to strings is expensive and increases file size several fold. In scientific computing, where the data is often only a few types of very large arrays, raw data is not that hard to convert into forms that enable data processing, and short descriptors can be written once to aid the conversion. However when you have more complex data, like when writing a time series with a variety of types, this shifts a lot of effort to the data analysis side. It also makes software updates dangerous because of the required consistency between logging and analysis.

This has come up several times for me when working on embedded projects, and I don’t have a great solution but I can share what I’ve settled on in C/C++. Please comment with suggestions if you have ‘killed’ this problem.

I have a microcontroller wired up to a number of sensors running at rates between 500Hz and 1Hz and I’m logging the information to a SD card. Writes are batched with a double buffering system such that they occur less often than the worst case SD card latency (~250ms). Each buffer has a series of structs prefixed by a typecode.

The first pain point here is the typecode. Given a list of packet definitions, how do you enforce that all the typecodes are unique, so that the most impaired contributor won’t mess anything up?

With python, you can do something like the following:

import inspect

class packet_A:

    def __init__(a):
        self.a = a

class packet_B:

    def __init__(a, b):
        self.a = a
        self.b = b

# After	the definitions, register typecodes                                                                                                                                                                          
def filter(name):
    return "packet" in name

class_list = [obj for obj in globals().values() if
              inspect.isclass(obj) and
              obj.__module__ == __name__ and
              filter(obj.__name__)]

counter = 1
for cls in class_list:
    cls.typecode = counter
    counter += 1

When anyone adds a packet definition the packet will get its own unique typecode. Modifying class properties at runtime can make debugging harder, but a search for ‘typecode’ will reveal this file and hopefully make it clear what is going wrong.

Until reflection comes to C++, this style of thing is off the table. One can do something similar with macros and the compiler extension __counter__, but it’s unclear enough that doing things by hand is better. I just make sure that the pattern is visible in header file where packets are defined:

struct packet_rtc {
  struct timer_time t;
  RTC_TimeTypeDef sTime;
  RTC_DateTypeDef sDate;
};

struct packet_vbatt {
  struct timer_time t;
  uint16_t vbatt_cnts;
};

struct packet_imu {
  struct timer_time t;
  uint16_t readings_cnts[4];
};

__inline__ uint8_t typecode(struct packet_rtc p) {
  return 1;
};

__inline__ uint8_t typecode(struct packet_vbatt) {
  return 2;
}

__inline__ uint8_t typecode(struct packet_imu) {
  return 3;
}

Now that that is ‘solved’, is there an automated way to find the structure of these packets (member offsets and sizes)? I’ve been fighting with this for a bit and playing with different C++ reflection libraries; ideally all the work would be done at compile time, as the information is all available then. In the end I settled on automatically generating the descriptor generating file using a python c++ header parsing library:

import cxxheaderparser.simple

header_text = open("Core/Inc/packet.hpp","r").read().replace("__inline__", "")

data = cxxheaderparser.simple.parse_string(header_text)

work_string = ""

for class_type in data.namespace.classes:
    class_name = class_type.class_decl.typename.segments[0].name
    to_format_string = "%s %d"
    format_args = f'"{class_name}", sizeof({class_name})'
    for field in class_type.fields:
        if isinstance(field.type, cxxheaderparser.types.Array):
            element_type = field.type.array_of.typename.segments[0].name
            element_num = field.type.size.tokens[0].value
            element_name = field.name
            to_format_string += " %s %d"
            format_args += f', "{element_type} {element_name}[{element_num}]", offsetof({class_name},{element_name})'
        elif isinstance(field.type, cxxheaderparser.types.Type):
            element_type = field.type.typename.segments[0].name
            element_name = field.name
            to_format_string += " %s %d"
            format_args += f', "{element_type} {element_name}", offsetof({class_name},{element_name})'
    to_format_string += "\\n\\r"
    work_string += f'    cx = snprintf(buff, sizeof(buff), "{to_format_string}", {format_args});\n'
    work_string += f'    CDC_Transmit_FS((uint8_t*)buff, cx);\n\n'

h_file_text = """                                                                                                                                                                                                    
void data_description();                                                                                                                                                                                             
"""

c_file_text = ("""#include "datadescriptor.hpp"                                                                                                                                                                      
#include <cstddef>                                                                                                                                                                                                   
#include <cstdio>                                                                                                                                                                                                    
#include "usbd_cdc_if.h"                                                                                                                                                                                             
#include "packet.hpp"                                                                                                                                                                                                
                                                                                                                                                                                                                     
void data_description() {                                                                                                                                                                                            
    char buff[200];                                                                                                                                                                                                  
    int cx;                                                                                                                                                                                                          
                                                                                                                                                                                                                     
""" + work_string +
"}\n")

with open("Core/Src/datadescriptor.cpp", "w") as f:
    f.write(c_file_text)

with open("Core/Inc/datadescriptor.hpp", "w") as f:
    f.write(h_file_text)

Then the function data_descriptor can be called on the embedded system to generate a string description with member types, names and offsets, which can be written at the beginning of the log and parsed to aid in data analysis.

I don’t find this a very satisfying solution, but it is better than writing out each string by hand for every packet type.

Light Propagation through Flesh with Blender Scripting

This blog post documents my attempts to set up a Blender pipeline to simulate the path of light through flesh, for studying PPG. I ultimately found out that Blender’s cycles rendering engine makes some unacceptable approximations with lighting.

Why

I’ve been back to working on biomonitoring projects (see https://akmemorph.us/biofeedback-smartwatch-update/, from 2019). Wearable SPO2 sensors/PPGs have become ubiquitous since then, and there have been some improvements in the field. There is research on where to put the photodiode compared to the light source (an advancement has been to use multiple photodiodes), and on different wavelengths of light, but there is also a conspicous absense of research on optics.

It is conceivable to me that controlling the distribution of radiance and the angular sensitivity of the photodiode would give better results. And given that phone cameras work at the size that they are, simple optics are probably possible in a watch formfactor. I’m also working on a wearable multi-sensor chest strap which will be less space constrained than a watch.

From looking at a variety of wearable PPG sensors, it appears like optics are very rare, though my Samsung Fit 3 appears to have a fresnel lens on it.

Plenty of research uses ray tracing to investigate PPG, without considering illumination or sensor optics (e.g. Optical Ray Tracing Simulation by Using Monte Carlo Method for
Reflectance-based Photoplethysmography Sensor in Human Skin and
Fingertip Model
). For this reason I wanted to get a ray tracing simulation set up.

I started using Blender for this project because it seemed like it has all the relevant physics (volumetric, anisotropic scattering and absorption, refraction and fresnel reflections). People have even modeled cameras in Blender for photorealistic lens abberations. Blender scripting is powerful and reasonable friendly, and I’ve used it professionally for digital twins and synthetic datasets for machine learning.

How

The following python script has a table of flesh optical properties from a paper, and adds a series of rectangular layers to an empty blender scene. It sets up a camera and a stand in light source. For multi-spectral work, the script could loop through optical properties at any given wavelength of interest.

# invoke me with	                                                                                                                                                                                                                                                                  
# ~> blender -P script.py
# Add -b to run in the background                                                                                                                                                                                                                                                       

import bpy
import math

class OpticalProperties:

    def __init__(self, IoR, scatter_pmm, absorption_pmm, anisotropy):
        self.IoR = IoR
        self.scatter_pmm = scatter_pmm
        self.absorption_pmm = absorption_pmm
	self.anisotropy = anisotropy

class Region:

    def __init__(self, name, thickness_mm, properties):
	self.name = name
        self.thickness_mm = thickness_mm
        self.properties = properties

light_distance_mm = 1.0

# From Optical Ray Tracing Simulation by Using Monte Carlo Method for ... computational and experiemntal research... 2022                                                                                                                                                         
# Uses Henyey-Greenstein for anisotropy                                                                                                                                                                                                                                           
regions = [Region('Air', light_distance_mm,
                  OpticalProperties(1., 0., 0., 0.)),
           Region('Epidermis', 0.2,
                  OpticalProperties(1.433, 10.7, 0.27, 0.79)),
           Region('Dermis_First_Layer', 0.2,
                  OpticalProperties(1.396, 18.7, 0.27, 0.81)),
           Region('Dermis_Plexus_Superficial', 0.2,
                  OpticalProperties(1.396, 19.2, 0.27, 0.82)),
           Region('Dermis_Second_Layer', 0.9,
                  OpticalProperties(1.396, 18.7, 0.27, 0.81)),
           Region('Plexus_Profundus', 0.6,
                  OpticalProperties(1.4, 22.5, 0.27, 0.78)),
           Region('Subcutaneous', 8.0,
                  OpticalProperties(1.37, 5.0, 0.003, 0.75))]

# Reset scene                                                                                                                                                                                                                                                                     
bpy.ops.wm.read_factory_settings(use_empty=True)

# Add flesh layers                                                                                                                                                                                                                                                                
starting_z_mm = 0
for i, region in enumerate(regions):
    bpy.ops.mesh.primitive_cube_add(size=1, location=(0, 0, starting_z_mm - 0.5 * region.thickness_mm), scale = (16, 16, region.thickness_mm))
    layer = bpy.context.view_layer.objects.active
    layer.name = "{:03d}_".format(i) + region.name
    mat = bpy.data.materials.new(name=region.name + "_mat")

    nodes = mat.node_tree.nodes
    links = mat.node_tree.links

    nodes.clear()

    volume_coefficients = nodes.new(type="ShaderNodeVolumeCoefficients")
    glass_bsdf = nodes.new(type="ShaderNodeBsdfGlass")
    output = nodes.new(type='ShaderNodeOutputMaterial')

    volume_coefficients.location = (0, -250)
    output.location = (300, 0)

    absorption = region.properties.absorption_pmm
    scatter = region.properties.scatter_pmm
    volume_coefficients.inputs['Absorption Coefficients'].default_value = (absorption, absorption, absorption)
    volume_coefficients.inputs['Scatter Coefficients'].default_value = (scatter, scatter, scatter)
    volume_coefficients.inputs['Anisotropy'].default_value = region.properties.anisotropy

    links.new(volume_coefficients.outputs['Volume'], output.inputs['Volume'])

    glass_bsdf.inputs['IOR'].default_value = region.properties.IoR

    links.new(glass_bsdf.outputs['BSDF'], output.inputs['Surface'])

    layer.data.materials.append(mat)

    starting_z_mm -= region.thickness_mm

# Add camera                                                                                                                                                                                                                                                                      
bpy.ops.object.camera_add(location=(0, 0, 15))
bpy.context.scene.camera = camera = bpy.context.view_layer.objects.active

# Add light                                                                                                                                                                                                                                                                       
bpy.ops.object.light_add(type='AREA', location=(3, 0, 0))
light = bpy.context.view_layer.objects.active
light.data.shape = 'RECTANGLE'
light.data.energy = 1000
light.data.size = 2.0
light.data.size_y = 1.0
light.data.spread = math.radians(60)

# Render settings                                                                                                                                                                                                                                                                 
scene = bpy.context.scene
scene.render.engine = 'CYCLES'
scene.render.filepath = "/tmp/output.png"
scene.render.resolution_x = 512
scene.render.resolution_y = 512

# Render                                                                                                                                                                                                                                                                          
#bpy.ops.render.render(write_still=True)

It sets up a shader block that looks like:

However, I noticed something suspicious. The image is different when the glass shader is missing from when it is there but with an IoR (Index of Refraction) set to 1. This raised the question: how is Blender doing the IoR anyway? A 2D plane with an refractive shader was visible and changed the image viewed through it, which is wrong, but it is conceivable to me that one could have rays remember the IoR of the medium they are in and update when crossing surfaces in a way that would be physically correct but wouldn’t work for non-manifold objects.

A pure refractive cylinder seems to have a plausible shadow and a plausible lensing effect. However, with an IoR of 1, the cylinder disappears (as it should) but leaves a shadow (as it shouldn’t).

I did more research into this, and it appears like it’s just a problem of the approximations Cycle’s path tracing make regarding lighting (see https://docs.blender.org/manual/nb/4.0/render/cycles/optimizations/reducing_noise.html#render-cycles-reducing-noise-glass-and-transp-shadows). There appear to be good hacks around it to improve the appearance in rendering edge cases like this, but no hacks that will give the physically realistic behavior that I want.

I’m moving to a real ray tracer (http://mcx.space/ or https://mitsuba.readthedocs.io/en/stable/index.html), but that’s a future post.