Combine two or more DensityPlot3Ds with different color schemes

I have two DensityPlot3D objects where the opacity and coloring are both physically meaningful, and I'm looking for a way to overlap them so that the resulting opacity/color is visually combined in an intuitive way. Thus, I'm pretty sure that combining them into one DensityPlot3D is not an option.

For example, consider the eigenmodes in the following situation:

reg = DiscretizeRegion[
ImplicitRegion[-1 <= x <= 1 && -1 <= y <= 1 && -1 <= z <= 1, {x, y, z}]]

{vals, funcs} =
NDEigensystem[
{Laplacian[u[x, y, z], {x, y, z}], DirichletCondition[u[x, y, z] == 0, True]},
u, {x, y, z} ∈ reg, 4]

g1 = DensityPlot3D[funcs[[2]][x, y, z], {x, y, z} ∈ reg,
ColorFunction -> "SunsetColors", Boxed -> False, Axes -> False, ViewPoint -> {3, 2, 3}];

g2 = DensityPlot3D[funcs[[3]][x, y, z], {x, y, z} ∈ reg,
ColorFunction -> "RoseColors", Boxed -> False, Axes -> False, ViewPoint -> {3, 2, 3}];


What I want is to be able to combine these two graphics in such a way that in the regions where the densities overlap the combined color/opacity is somehow the sum or product or some other meaningful combination of the individual densities, but maintaining the different color schemes for the two modes (ie, not one unified DensityPlot3D with a common color function and opacity scaling).

I could see calculating an effective color and opacity value at discrete points in the region, but there does not seem to be a way to define a color function or opacity function that depend on position as well as "value".

EDIT: Things that don't work:

DensityPlot3D[funcs[[2]][x, y, z], {x, y, z} ∈ reg, ColorFunction -> Function[{x, y, z, f}, <anything>]]


(fails because ColorFunction doesn't work that way for DensityPlot3D)

DensityPlot3D[{funcs[[2]][x, y, z],funcs[[3]][x, y, z]}, {x, y, z} ∈ reg]


(displays nothing)

Show[g1, g2]


(only shows g1, because it rasterizes first)

Blend[{g1,g2}]


(rasterizes and then applies the same opacity at every pixel)

• g1 and g2 appear identical to me in your code. Did you mean to select a different part of funcs for one of them? Also, since your region is a simple axis-aligned cube, is it faster to use reg = Cuboid[{-1, -1, -1}, {1, 1, 1}]? – MarcoB Apr 2 '19 at 22:51
• Can you not define a ColorFunction using both position and density values? ColorFunction -> Function[{x, y, z, f}, yourfunc[x,y,z,f]] – MelaGo Apr 2 '19 at 23:38
• Edited to correct g2. As for ColorFunction, I get all sorts of errors when trying to do what you suggest in DensityPlot3D. Maybe there's another way to do it with ListPlot3D? I'm not smart enough, though. – djphd Apr 3 '19 at 1:25
• @djphd What are "all sorts of errors"? Could you add a few examples of what you tried together with the errors to the question? – Lukas Lang Apr 8 '19 at 11:44
• I would try to construct an Image3D directly instead of using DensityPlot3D. Then I could produce RGB+alpha quadruplets manually and control how the colour is created. I would no longer need to compute the colour based on a single value: I could compute it based on two values (which is what you are asking). – Szabolcs Apr 9 '19 at 15:05

Here's an answer. It's a long workaround, so I'll happily accept another answer that takes advantage of simple, built-in functionality. First I construct tables of (rescaled/normalized) function values for funcs[[2]] and funcs[[3]] called r1 and r2, respectively. Then I make my own table of RGBColor values with opacity as the sum of squares of the function values (there's probably a more elegant way to do that than what I came up with, too).
colorTable = Table[RGBColor[