diff --git a/stuff/config/current.txt b/stuff/config/current.txt
index 528f706..a2b8954 100644
--- a/stuff/config/current.txt
+++ b/stuff/config/current.txt
@@ -1174,6 +1174,27 @@
   <item>"STD_iwa_PNPerspectiveFx.alpha_rendering"	"Alpha Rendering"</item>
   <item>"STD_iwa_PNPerspectiveFx.waveHeight"		"Wave Height"</item>
   
+  <item>"STD_iwa_SoapBubbleFx"	"SoapBubble Iwa"	</item>
+  <item>"STD_iwa_SoapBubbleFx.intensity"	"Intensity"	</item>
+  <item>"STD_iwa_SoapBubbleFx.refractiveIndex"	"Refractive Index"	</item>
+  <item>"STD_iwa_SoapBubbleFx.thickMax"	"Thick Max"	</item>
+  <item>"STD_iwa_SoapBubbleFx.thickMin"	"Thick Min"	</item>
+  <item>"STD_iwa_SoapBubbleFx.RGamma"	"R Gamma"	</item>
+  <item>"STD_iwa_SoapBubbleFx.GGamma"	"G Gamma"	</item>
+  <item>"STD_iwa_SoapBubbleFx.BGamma"	"B Gamma"	</item>
+  <item>"STD_iwa_SoapBubbleFx.binarizeThresold"	"Threshold"	</item>
+  <item>"STD_iwa_SoapBubbleFx.shapeAspectRatio"	"Shape Aspect Ratio"	</item>
+  <item>"STD_iwa_SoapBubbleFx.blurRadius"	"Blur Radius"	</item>
+  <item>"STD_iwa_SoapBubbleFx.blurPower"	"Power"	</item>
+  <item>"STD_iwa_SoapBubbleFx.normalSampleDistance"	"Sample Distance"	</item>
+  <item>"STD_iwa_SoapBubbleFx.noiseSubDepth"	"Sub Depth"	</item>
+  <item>"STD_iwa_SoapBubbleFx.noiseResolutionS"	"S Resolution"	</item>
+  <item>"STD_iwa_SoapBubbleFx.noiseResolutionT"	"T Resolution"	</item>
+  <item>"STD_iwa_SoapBubbleFx.noiseSubCompositeRatio"	"Sub Amplitude Ratio"	</item>
+  <item>"STD_iwa_SoapBubbleFx.noiseEvolution"	"Evolution"	</item>
+  <item>"STD_iwa_SoapBubbleFx.noiseDepthMixRatio"	"Noise to Depth"	</item>
+  <item>"STD_iwa_SoapBubbleFx.noiseThicknessMixRatio"	"Noise to Thickness"	</item>
+
  <!------------------------------ Tiled Particles Iwa ------------------------------------------->
 
   <item>STD_iwa_TiledParticlesFx "Tiled Particles Iwa" </item>
diff --git a/stuff/profiles/layouts/fxs/STD_iwa_SoapBubbleFx.xml b/stuff/profiles/layouts/fxs/STD_iwa_SoapBubbleFx.xml
new file mode 100644
index 0000000..1180dbc
--- /dev/null
+++ b/stuff/profiles/layouts/fxs/STD_iwa_SoapBubbleFx.xml
@@ -0,0 +1,33 @@
+<fxlayout>
+  <page name="Color and Shape">
+    <vbox>
+      <separator label="Bubble Color"/>
+      <control>intensity</control>
+      <control>refractiveIndex</control>
+      <control>thickMax</control>
+      <control>thickMin</control>
+      <control>RGamma</control>
+      <control>GGamma</control>
+      <control>BGamma</control>
+    </vbox>
+    <vbox>
+      <separator label="Shape"/>
+      <control>binarizeThresold</control>
+      <control>shpeAspectRatio</control>
+      <control>blurRadius</control>
+      <control>blurPower</control>
+    </vbox>
+  </page>
+  <page name="Noise">
+    <vbox>
+      <control>normalSampleDistance</control>
+      <control>noiseSubDepth</control>
+      <control>noiseResolutionS</control>
+      <control>noiseResolutionT</control>
+      <control>noiseSubCompositeRatio</control>
+      <control>noiseEvolution</control>
+      <control>noiseDepthMixRatio</control>
+      <control>noiseThicknessMixRatio</control>
+    </vbox>
+  </page>
+</fxlayout>
diff --git a/stuff/profiles/layouts/fxs/fxs.lst b/stuff/profiles/layouts/fxs/fxs.lst
index 509e17e..2ff0aba 100644
--- a/stuff/profiles/layouts/fxs/fxs.lst
+++ b/stuff/profiles/layouts/fxs/fxs.lst
@@ -115,8 +115,9 @@
     STD_inoFogFx
     STD_glowFx
     STD_lightSpotFx
-    STD_raylitFx 
+    STD_raylitFx
     STD_iwa_SpectrumFx
+    STD_iwa_SoapBubbleFx
     STD_targetSpotFx
   </Light>
   <Matte>
diff --git a/toonz/sources/stdfx/CMakeLists.txt b/toonz/sources/stdfx/CMakeLists.txt
index 3e45b43..3f78804 100644
--- a/toonz/sources/stdfx/CMakeLists.txt
+++ b/toonz/sources/stdfx/CMakeLists.txt
@@ -71,6 +71,7 @@ set(HEADERS
     iwa_noise1234.h
     iwa_fresnel.h
     iwa_pnperspectivefx.h
+    iwa_soapbubblefx.h
 )
 
 set(SOURCES
@@ -244,6 +245,7 @@ set(SOURCES
     iwa_simplexnoise.cpp
     iwa_noise1234.cpp
     iwa_pnperspectivefx.cpp
+    iwa_soapbubblefx.cpp
 )
 
 set(OBJCSOURCES
diff --git a/toonz/sources/stdfx/iwa_soapbubblefx.cpp b/toonz/sources/stdfx/iwa_soapbubblefx.cpp
new file mode 100644
index 0000000..24249f6
--- /dev/null
+++ b/toonz/sources/stdfx/iwa_soapbubblefx.cpp
@@ -0,0 +1,714 @@
+/*------------------------------------
+Iwa_SoapBubbleFx
+Generates thin film interference colors from two reference images;
+one is for thickness and the other one is for shape or normal vector
+distribution of the film.
+Inherits Iwa_SpectrumFx.
+------------------------------------*/
+
+#include "iwa_soapbubblefx.h"
+#include "iwa_cie_d65.h"
+#include "iwa_xyz.h"
+
+#include <QList>
+#include <QPoint>
+#include <QSize>
+
+namespace {
+const float PI = 3.14159265f;
+
+#define INF 1e20 /* less than FLT_MAX */
+
+/* dt of 1d function using squared distance */
+static float* dt(float* f, int n, float a = 1.0f) {
+  float* d = new float[n];
+  int* v   = new int[n];
+  float* z = new float[n + 1];
+  /* index of rightmost parabola in lower envelope */
+  int k = 0;
+  /* locations of parabolas in lower envelope */
+  v[0] = 0;
+  /* locations of boundaries between parabolas */
+  z[0] = -INF;
+  z[1] = +INF;
+  /* compute lower envelope */
+  for (int q = 1; q <= n - 1; q++) {
+    /* compute intersection */
+    float s =
+        ((f[q] / a + q * q) - (f[v[k]] / a + v[k] * v[k])) / (2 * q - 2 * v[k]);
+    while (s <= z[k]) {
+      k--;
+      s = ((f[q] / a + q * q) - (f[v[k]] / a + v[k] * v[k])) /
+          (2 * q - 2 * v[k]);
+    }
+    k++;
+    v[k]     = q;
+    z[k]     = s;
+    z[k + 1] = +INF;
+  }
+  k = 0;
+  /* fill in values of distance transform */
+  for (int q = 0; q <= n - 1; q++) {
+    while (z[k + 1] < q) k++;
+    d[q] = a * (q - v[k]) * (q - v[k]) + f[v[k]];
+  }
+
+  delete[] v;
+  delete[] z;
+  return d;
+}
+}
+
+//------------------------------------
+
+Iwa_SoapBubbleFx::Iwa_SoapBubbleFx()
+    : Iwa_SpectrumFx()
+    , m_binarize_threshold(0.5)
+    , m_shape_aspect_ratio(1.0)
+    , m_blur_radius(5.0)
+    , m_blur_power(0.5)
+    , m_normal_sample_distance(1)
+    , m_noise_sub_depth(3)
+    , m_noise_resolution_s(18.0)
+    , m_noise_resolution_t(5.0)
+    , m_noise_sub_composite_ratio(0.5)
+    , m_noise_evolution(0.0)
+    , m_noise_depth_mix_ratio(0.05)
+    , m_noise_thickness_mix_ratio(0.05) {
+  removeInputPort("Source");
+  removeInputPort("Light"); /* not used */
+  addInputPort("Thickness", m_input);
+  addInputPort("Shape", m_shape);
+  addInputPort("Depth", m_depth);
+
+  bindParam(this, "binarizeThresold", m_binarize_threshold);
+  bindParam(this, "shapeAspectRatio", m_shape_aspect_ratio);
+  bindParam(this, "blurRadius", m_blur_radius);
+  bindParam(this, "blurPower", m_blur_power);
+  bindParam(this, "normalSampleDistance", m_normal_sample_distance);
+  bindParam(this, "noiseSubDepth", m_noise_sub_depth);
+  bindParam(this, "noiseResolutionS", m_noise_resolution_s);
+  bindParam(this, "noiseResolutionT", m_noise_resolution_t);
+  bindParam(this, "noiseSubCompositeRatio", m_noise_sub_composite_ratio);
+  bindParam(this, "noiseEvolution", m_noise_evolution);
+  bindParam(this, "noiseDepthMixRatio", m_noise_depth_mix_ratio);
+  bindParam(this, "noiseThicknessMixRatio", m_noise_thickness_mix_ratio);
+
+  m_binarize_threshold->setValueRange(0.01, 0.99);
+  m_shape_aspect_ratio->setValueRange(0.2, 5.0);
+  m_blur_radius->setMeasureName("fxLength");
+  m_blur_radius->setValueRange(0.0, 25.0);
+  m_blur_power->setValueRange(0.01, 5.0);
+
+  m_normal_sample_distance->setValueRange(1, 20);
+  m_noise_sub_depth->setValueRange(1, 5);
+  m_noise_resolution_s->setValueRange(1.0, 40.0);
+  m_noise_resolution_t->setValueRange(1.0, 20.0);
+  m_noise_sub_composite_ratio->setValueRange(0.0, 5.0);
+  m_noise_depth_mix_ratio->setValueRange(0.0, 1.0);
+  m_noise_thickness_mix_ratio->setValueRange(0.0, 1.0);
+}
+
+//------------------------------------
+
+void Iwa_SoapBubbleFx::doCompute(TTile& tile, double frame,
+                                 const TRenderSettings& settings) {
+  if (!m_input.isConnected()) return;
+  if (!m_shape.isConnected() && !m_depth.isConnected()) return;
+
+  TDimensionI dim(tile.getRaster()->getLx(), tile.getRaster()->getLy());
+  TRectD bBox(tile.m_pos, TPointD(dim.lx, dim.ly));
+
+  /* soap bubble color map */
+  TRasterGR8P bubbleColor_ras(sizeof(float3) * 256 * 256, 1);
+  bubbleColor_ras->lock();
+  float3* bubbleColor_p = (float3*)bubbleColor_ras->getRawData();
+  calcBubbleMap(bubbleColor_p, frame, true);
+
+  /* depth map */
+  TRasterGR8P depth_map_ras(sizeof(float) * dim.lx * dim.ly, 1);
+  depth_map_ras->lock();
+  float* depth_map_p = (float*)depth_map_ras->getRawData();
+
+  /* if the depth image is connected, use it */
+  if (m_depth.isConnected()) {
+    TTile depth_tile;
+    m_depth->allocateAndCompute(depth_tile, bBox.getP00(), dim,
+                                tile.getRaster(), frame, settings);
+    TRasterP depthRas = depth_tile.getRaster();
+    depthRas->lock();
+
+    TRaster32P depthRas32 = (TRaster32P)depthRas;
+    TRaster64P depthRas64 = (TRaster64P)depthRas;
+    {
+      if (depthRas32)
+        convertToBrightness<TRaster32P, TPixel32>(depthRas32, depth_map_p, dim);
+      else if (depthRas64)
+        convertToBrightness<TRaster64P, TPixel64>(depthRas64, depth_map_p, dim);
+    }
+    depthRas->unlock();
+  }
+  /* or, use the shape image to obtain pseudo depth */
+  else { /* m_shape.isConnected */
+    /* obtain shape image */
+    TTile shape_tile;
+    {
+      TRaster32P tmp(1, 1);
+      m_shape->allocateAndCompute(shape_tile, bBox.getP00(), dim, tmp, frame,
+                                  settings);
+    }
+    processShape(frame, shape_tile, depth_map_p, dim, settings);
+  }
+
+  /* compute the thickness input and temporarily store to the tile */
+  m_input->compute(tile, frame, settings);
+
+  TRasterGR8P thickness_map_ras(sizeof(float) * dim.lx * dim.ly, 1);
+  thickness_map_ras->lock();
+  float* thickness_map_p = (float*)thickness_map_ras->getRawData();
+  TRasterP thicknessRas  = tile.getRaster();
+  TRaster32P ras32       = (TRaster32P)thicknessRas;
+  TRaster64P ras64       = (TRaster64P)thicknessRas;
+  {
+    if (ras32)
+      convertToBrightness<TRaster32P, TPixel32>(ras32, thickness_map_p, dim);
+    else if (ras64)
+      convertToBrightness<TRaster64P, TPixel64>(ras64, thickness_map_p, dim);
+  }
+
+  /* process noise */
+  processNoise(thickness_map_p, depth_map_p, dim, frame, settings);
+
+  if (ras32)
+    convertToRaster<TRaster32P, TPixel32>(ras32, thickness_map_p, depth_map_p,
+                                          dim, bubbleColor_p);
+  else if (ras64)
+    convertToRaster<TRaster64P, TPixel64>(ras64, thickness_map_p, depth_map_p,
+                                          dim, bubbleColor_p);
+
+  thickness_map_ras->unlock();
+  depth_map_ras->unlock();
+  bubbleColor_ras->unlock();
+}
+
+//------------------------------------
+
+template <typename RASTER, typename PIXEL>
+void Iwa_SoapBubbleFx::convertToBrightness(const RASTER srcRas, float* dst,
+                                           TDimensionI dim) {
+  float* dst_p = dst;
+  for (int j = 0; j < dim.ly; j++) {
+    PIXEL* pix = srcRas->pixels(j);
+    for (int i = 0; i < dim.lx; i++, dst_p++, pix++) {
+      float r = (float)pix->r / (float)PIXEL::maxChannelValue;
+      float g = (float)pix->g / (float)PIXEL::maxChannelValue;
+      float b = (float)pix->b / (float)PIXEL::maxChannelValue;
+      /* brightness */
+      *dst_p = 0.298912f * r + 0.586611f * g + 0.114478f * b;
+    }
+  }
+}
+
+//------------------------------------
+
+template <typename RASTER, typename PIXEL>
+void Iwa_SoapBubbleFx::convertToRaster(const RASTER ras, float* thickness_map_p,
+                                       float* depth_map_p, TDimensionI dim,
+                                       float3* bubbleColor_p) {
+  float* depth_p     = depth_map_p;
+  float* thickness_p = thickness_map_p;
+  for (int j = 0; j < dim.ly; j++) {
+    PIXEL* pix = ras->pixels(j);
+    for (int i = 0; i < dim.lx; i++, depth_p++, thickness_p++, pix++) {
+      float alpha = (float)pix->m / PIXEL::maxChannelValue;
+      if (alpha == 0.0f) /* no change for the transparent pixels */
+        continue;
+
+      float coordinate[2];
+      coordinate[0] = 256.0f * std::min(1.0f, *depth_p);
+      coordinate[1] = 256.0f * std::min(1.0f, *thickness_p);
+
+      int neighbors[2][2];
+
+      /* interpolate sampling */
+      if (coordinate[0] <= 0.5f)
+        neighbors[0][0] = 0;
+      else
+        neighbors[0][0] = (int)std::floor(coordinate[0] - 0.5f);
+      if (coordinate[0] >= 255.5f)
+        neighbors[0][1] = 255;
+      else
+        neighbors[0][1] = (int)std::floor(coordinate[0] + 0.5f);
+      if (coordinate[1] <= 0.5f)
+        neighbors[1][0] = 0;
+      else
+        neighbors[1][0] = (int)std::floor(coordinate[1] - 0.5f);
+      if (coordinate[1] >= 255.5f)
+        neighbors[1][1] = 255;
+      else
+        neighbors[1][1] = (int)std::floor(coordinate[1] + 0.5f);
+
+      float interp_ratio[2];
+      interp_ratio[0] = coordinate[0] - 0.5f - std::floor(coordinate[0] - 0.5f);
+      interp_ratio[1] = coordinate[1] - 0.5f - std::floor(coordinate[1] - 0.5f);
+
+      float3 nColors[4] = {
+          bubbleColor_p[neighbors[0][0] * 256 + neighbors[1][0]],
+          bubbleColor_p[neighbors[0][1] * 256 + neighbors[1][0]],
+          bubbleColor_p[neighbors[0][0] * 256 + neighbors[1][1]],
+          bubbleColor_p[neighbors[0][1] * 256 + neighbors[1][1]]};
+
+      float3 color =
+          nColors[0] * (1.0f - interp_ratio[0]) * (1.0f - interp_ratio[1]) +
+          nColors[1] * interp_ratio[0] * (1.0f - interp_ratio[1]) +
+          nColors[2] * (1.0f - interp_ratio[0]) * interp_ratio[1] +
+          nColors[3] * interp_ratio[0] * interp_ratio[1];
+
+      /* clamp */
+      float val = color.x * (float)PIXEL::maxChannelValue + 0.5f;
+      pix->r = (typename PIXEL::Channel)((val > (float)PIXEL::maxChannelValue)
+                                             ? (float)PIXEL::maxChannelValue
+                                             : val);
+      val    = color.y * (float)PIXEL::maxChannelValue + 0.5f;
+      pix->g = (typename PIXEL::Channel)((val > (float)PIXEL::maxChannelValue)
+                                             ? (float)PIXEL::maxChannelValue
+                                             : val);
+      val    = color.z * (float)PIXEL::maxChannelValue + 0.5f;
+      pix->b = (typename PIXEL::Channel)((val > (float)PIXEL::maxChannelValue)
+                                             ? (float)PIXEL::maxChannelValue
+                                             : val);
+    }
+  }
+}
+
+//------------------------------------
+
+void Iwa_SoapBubbleFx::processShape(double frame, TTile& shape_tile,
+                                    float* depth_map_p, TDimensionI dim,
+                                    const TRenderSettings& settings) {
+  TRaster32P shapeRas = shape_tile.getRaster();
+  shapeRas->lock();
+
+  /* binarize the shape image */
+  TRasterGR8P binarized_ras(sizeof(char) * dim.lx * dim.ly, 1);
+  binarized_ras->lock();
+  char* binarized_p = (char*)binarized_ras->getRawData();
+
+  TRasterGR8P distance_ras(sizeof(float) * dim.lx * dim.ly, 1);
+  distance_ras->lock();
+  float* distance_p = (float*)distance_ras->getRawData();
+
+  float binarize_thres = (float)m_binarize_threshold->getValue(frame);
+
+  do_binarize(shapeRas, binarized_p, binarize_thres, distance_p, dim);
+
+  shapeRas->unlock();
+
+  do_distance_transform(distance_p, dim, frame);
+
+  /* create blur filter */
+  float blur_radius = (float)m_blur_radius->getValue(frame) *
+                      std::sqrt(std::abs((float)settings.m_affine.det()));
+
+  /* if blur radius is 0, set the distance image to the depth image as-is */
+  if (blur_radius == 0.0f) {
+    float power      = (float)m_blur_power->getValue(frame);
+    float* tmp_depth = depth_map_p;
+    float* tmp_dist  = distance_p;
+    char* bin_p      = binarized_p;
+    for (int i = 0; i < dim.lx * dim.ly;
+         i++, tmp_depth++, tmp_dist++, bin_p++) {
+      if (*bin_p == 0)
+        *tmp_depth = 0.0f;
+      else
+        *tmp_depth = 1.0f - std::pow(*tmp_dist, power);
+    }
+    distance_ras->unlock();
+    binarized_ras->unlock();
+    return;
+  }
+
+  int blur_filter_size = (int)std::floor(blur_radius) * 2 + 1;
+  TRasterGR8P blur_filter_ras(
+      sizeof(float) * blur_filter_size * blur_filter_size, 1);
+  blur_filter_ras->lock();
+  float* blur_filter_p = (float*)blur_filter_ras->getRawData();
+
+  do_createBlurFilter(blur_filter_p, blur_filter_size, blur_radius);
+
+  /* blur filtering, normarize & power */
+  do_applyFilter(depth_map_p, dim, distance_p, binarized_p, blur_filter_p,
+                 blur_filter_size, frame);
+
+  distance_ras->unlock();
+  binarized_ras->unlock();
+  blur_filter_ras->unlock();
+}
+
+//------------------------------------
+
+void Iwa_SoapBubbleFx::do_binarize(TRaster32P srcRas, char* dst_p, float thres,
+                                   float* distance_p, TDimensionI dim) {
+  TPixel32::Channel channelThres =
+      (TPixel32::Channel)(thres * (float)TPixel32::maxChannelValue);
+  char* tmp_p     = dst_p;
+  float* tmp_dist = distance_p;
+  for (int j = 0; j < dim.ly; j++) {
+    TPixel32* pix = srcRas->pixels(j);
+    for (int i = 0; i < dim.lx; i++, pix++, tmp_p++, tmp_dist++) {
+      (*tmp_p)    = (pix->m > channelThres) ? 1 : 0;
+      (*tmp_dist) = (*tmp_p == 1) ? INF : 0.0f;
+    }
+  }
+}
+
+//------------------------------------
+
+void Iwa_SoapBubbleFx::do_createBlurFilter(float* dst_p, int size,
+                                           float radius) {
+  float radius2 = radius * radius;
+  float* tmp_p  = dst_p;
+  float sum     = 0.0f;
+  int rad       = (size - 1) / 2;
+  for (int j = -rad; j <= rad; j++) {
+    for (int i = -rad; i <= rad; i++, tmp_p++) {
+      float length2 = (float)i * (float)i + (float)j * (float)j;
+      /* out of range */
+      if (length2 >= radius2)
+        *tmp_p = 0.0f;
+      else {
+        /* normalize distace from the filter center, to 0-1 */
+        *tmp_p = 1.0f - std::sqrt(length2) / radius;
+        sum += *tmp_p;
+      }
+    }
+  }
+  /* normalize */
+  tmp_p = dst_p;
+  for (int i = 0; i < size * size; i++, tmp_p++) {
+    *tmp_p /= sum;
+  }
+}
+//------------------------------------
+
+void Iwa_SoapBubbleFx::do_applyFilter(float* depth_map_p, TDimensionI dim,
+                                      float* distance_p, char* binarized_p,
+                                      float* blur_filter_p,
+                                      int blur_filter_size, double frame) {
+  float power = (float)m_blur_power->getValue(frame);
+
+  memset(depth_map_p, 0, sizeof(float) * dim.lx * dim.ly);
+
+  int fil_margin = (blur_filter_size - 1) / 2;
+  float* dst_p   = depth_map_p;
+  char* bin_p    = binarized_p;
+  for (int j = 0; j < dim.ly; j++) {
+    for (int i = 0; i < dim.lx; i++, dst_p++, bin_p++) {
+      if (*bin_p == 0) continue;
+
+      float* fil_p = blur_filter_p;
+      for (int fy = j - fil_margin; fy <= j + fil_margin; fy++) {
+        if (fy < 0 || fy >= dim.ly) {
+          fil_p += blur_filter_size;
+          continue;
+        }
+        for (int fx = i - fil_margin; fx <= i + fil_margin; fx++, fil_p++) {
+          if (fx < 0 || fx >= dim.lx) continue;
+
+          *dst_p += *fil_p * distance_p[fy * dim.lx + fx];
+        }
+      }
+      /* power the value */
+      *dst_p = 1.0f - std::pow(*dst_p, power);
+    }
+  }
+}
+
+//------------------------------------
+
+void Iwa_SoapBubbleFx::processNoise(float* thickness_map_p, float* depth_map_p,
+                                    TDimensionI dim, double frame,
+                                    const TRenderSettings& settings) {
+  float noise_depth_mix_ratio = (float)m_noise_depth_mix_ratio->getValue(frame);
+  float noise_thickness_mix_ratio =
+      (float)m_noise_thickness_mix_ratio->getValue(frame);
+
+  /* If the noise ratio is 0, do nothing and return */
+  if (noise_depth_mix_ratio == 0.0f && noise_thickness_mix_ratio == 0.0f)
+    return;
+
+  int noise_sub_depth    = m_noise_sub_depth->getValue();
+  int noise_resolution_s = (int)m_noise_resolution_s->getValue(frame);
+  int noise_resolution_t = (int)m_noise_resolution_t->getValue(frame);
+  float noise_composite_ratio =
+      (float)m_noise_sub_composite_ratio->getValue(frame);
+  float noise_evolution = (float)m_noise_evolution->getValue(frame);
+
+  /* initialize the phase map */
+  QList<int> noise_amount;
+  QList<QSize> noise_base_resolution;
+  int whole_noise_amount = 0;
+
+  for (int layer = 0; layer < noise_sub_depth; layer++) {
+    /* noise resolution */
+    /* width: circumferential direction   height:distal direction */
+    QSize size;
+    size.setWidth(std::pow(2, layer) * noise_resolution_s);
+    size.setHeight(std::pow(2, layer) * noise_resolution_t + 1);
+    noise_base_resolution.append(size);
+    int amount = size.width() * size.height();
+    noise_amount.append(amount);
+    whole_noise_amount += amount;
+  }
+
+  float* noise_phases = new float[whole_noise_amount];
+  float* ph_p         = noise_phases;
+
+  srand(0);
+  /* Set the phase differences (0-2��) */
+  for (int i = 0; i < whole_noise_amount; i++, ph_p++) {
+    *ph_p = (float)rand() / (float)RAND_MAX * 2.0f * PI;
+  }
+
+  /* make noise base */
+  /* compute composite ratio of each layer */
+  QList<float> comp_ratios;
+  comp_ratios.append(10.0f);
+  float ratio_sum = 10.0f;
+  for (int i = 1; i < noise_sub_depth; i++) {
+    comp_ratios.append(comp_ratios.last() * noise_composite_ratio);
+    ratio_sum += comp_ratios.last();
+  }
+  /* normalize */
+  for (int i = 0; i < noise_sub_depth; i++) comp_ratios[i] /= ratio_sum;
+
+  float* noise_base = new float[whole_noise_amount];
+
+  float* nb_p = noise_base;
+  ph_p        = noise_phases;
+
+  /* for each sub-noise layer */
+  for (int layer = 0; layer < noise_sub_depth; layer++) {
+    float tmp_evolution = noise_evolution * (float)(layer + 1);
+    for (int i = 0; i < noise_amount[layer]; i++, nb_p++, ph_p++) {
+      *nb_p = comp_ratios[layer] * (cosf(tmp_evolution + *ph_p) / 2.0f + 0.5f);
+    }
+  }
+  delete[] noise_phases;
+
+  TRasterGR8P norm_angle_ras(sizeof(float) * dim.lx * dim.ly, 1);
+  norm_angle_ras->lock();
+  float* norm_angle_p = (float*)norm_angle_ras->getRawData();
+
+  calc_norm_angle(norm_angle_p, depth_map_p, dim, settings.m_shrinkX);
+
+  TRasterGR8P noise_map_ras(sizeof(float) * dim.lx * dim.ly, 1);
+  noise_map_ras->lock();
+  float* noise_map_p = (float*)noise_map_ras->getRawData();
+
+  make_noise_map(noise_map_p, depth_map_p, norm_angle_p, dim, noise_amount,
+                 noise_base_resolution, noise_sub_depth, noise_base);
+
+  norm_angle_ras->unlock();
+  delete[] noise_base;
+
+  /* composite with perlin noise */
+  add_noise(thickness_map_p, depth_map_p, dim, noise_map_p,
+            noise_thickness_mix_ratio, noise_depth_mix_ratio);
+
+  noise_map_ras->unlock();
+}
+
+//------------------------------------
+
+void Iwa_SoapBubbleFx::calc_norm_angle(float* norm_angle_p, float* depth_map_p,
+                                       TDimensionI dim, int shrink) {
+  struct Locals {
+    TDimensionI _dim;
+    const float* _depth_p;
+    float data(int x, int y) {
+      if (x < 0 || _dim.lx <= x || y < 0 || _dim.ly <= y) return 0.0f;
+      return _depth_p[y * _dim.lx + x];
+    }
+  } locals = {dim, depth_map_p};
+
+  int sampleDistance =
+      std::max(1, m_normal_sample_distance->getValue() / shrink);
+  float* dst_p = norm_angle_p;
+
+  for (int j = 0; j < dim.ly; j++) {
+    int sample_y[2]                  = {j - sampleDistance, j + sampleDistance};
+    if (sample_y[0] < 0) sample_y[0] = 0;
+    if (sample_y[1] >= dim.ly) sample_y[1] = dim.ly - 1;
+
+    for (int i = 0; i < dim.lx; i++, norm_angle_p++) {
+      int sample_x[2] = {i - sampleDistance, i + sampleDistance};
+      if (sample_x[1] >= dim.lx) sample_x[1] = dim.lx - 1;
+      if (sample_x[0] < 0) sample_x[0]       = 0;
+
+      float gradient[2];
+      gradient[0] =
+          (locals.data(sample_x[0], j) - locals.data(sample_x[1], j)) /
+          (float)(sample_x[0] - sample_x[1]);
+      gradient[1] =
+          (locals.data(i, sample_y[0]) - locals.data(i, sample_y[1])) /
+          (float)(sample_y[0] - sample_y[1]);
+
+      if (gradient[0] == 0.0f && gradient[1] == 0.0f)
+        *norm_angle_p = 0.0f;
+      else /* normalize value range to 0-1 */
+        *norm_angle_p =
+            0.5f + std::atan2(gradient[0], gradient[1]) / (2.0f * PI);
+    }
+  }
+}
+
+//------------------------------------
+
+void Iwa_SoapBubbleFx::make_noise_map(float* noise_map_p, float* depth_map_p,
+                                      float* norm_angle_p, TDimensionI dim,
+                                      const QList<int>& noise_amount,
+                                      const QList<QSize>& noise_base_resolution,
+                                      int noise_sub_depth, float* noise_base) {
+  float* dst_p   = noise_map_p;
+  float* depth_p = depth_map_p;
+  float* norm_p  = norm_angle_p;
+
+  for (int j = 0; j < dim.ly; j++) {
+    for (int i = 0; i < dim.lx; i++, dst_p++, depth_p++, norm_p++) {
+      /* Obtain coordinate */
+      /* circumferential direction */
+      float tmp_s = (*norm_p);
+      /* distal direction */
+      float tmp_t = std::min(1.0f, *depth_p);
+
+      /* accumulate noise values */
+      *dst_p                  = 0.0f;
+      float* noise_layer_base = noise_base;
+      for (int layer = 0; layer < noise_sub_depth; layer++) {
+        /* obtain pseudo polar coords */
+        QSize reso = noise_base_resolution.at(layer);
+        float polar_s =
+            tmp_s * (float)(reso.width()); /* because it is circumferential */
+        float polar_t = tmp_t * (float)(reso.height() - 1);
+
+        /* first, compute circumferential position and ratio */
+        int neighbor_s[2];
+        neighbor_s[0] = (int)std::floor(polar_s);
+        neighbor_s[1] = neighbor_s[0] + 1;
+        if (neighbor_s[0] == reso.width()) neighbor_s[0] = 0;
+        if (neighbor_s[1] >= reso.width()) neighbor_s[1] = 0;
+        float ratio_s = polar_s - std::floor(polar_s);
+
+        /* second, compute distal position and ratio */
+        int neighbor_t[2];
+        neighbor_t[0] = (int)std::floor(polar_t);
+        neighbor_t[1] = neighbor_t[0] + 1;
+        if (neighbor_t[1] == reso.height()) neighbor_t[1] -= 1;
+        float ratio_t = polar_t - std::floor(polar_t);
+
+        *dst_p += noise_interp(neighbor_s[0], neighbor_s[1], neighbor_t[0],
+                               neighbor_t[1], ratio_s, ratio_t,
+                               noise_layer_base, reso.width());
+
+        /* offset noise pointer */
+        noise_layer_base += noise_amount[layer];
+      }
+    }
+  }
+}
+
+//------------------------------------
+
+float Iwa_SoapBubbleFx::noise_interp(int left, int right, int bottom, int top,
+                                     float ratio_s, float ratio_t,
+                                     float* noise_layer_base, int noise_dim_x) {
+  struct Locals {
+    int _dim_x;
+    const float* _noise_p;
+    float data(int x, int y) { return _noise_p[y * _dim_x + x]; }
+  } locals = {noise_dim_x, noise_layer_base};
+
+  float c_ratio_s = (1.0f - cosf(ratio_s * PI)) * 0.5f;
+  float c_ratio_t = (1.0f - cosf(ratio_t * PI)) * 0.5f;
+
+  return locals.data(left, bottom) * (1.0f - c_ratio_s) * (1.0f - c_ratio_t) +
+         locals.data(right, bottom) * c_ratio_s * (1.0f - c_ratio_t) +
+         locals.data(left, top) * (1.0f - c_ratio_s) * c_ratio_t +
+         locals.data(right, top) * c_ratio_s * c_ratio_t;
+}
+
+//------------------------------------
+
+void Iwa_SoapBubbleFx::add_noise(float* thickness_map_p, float* depth_map_p,
+                                 TDimensionI dim, float* noise_map_p,
+                                 float noise_thickness_mix_ratio,
+                                 float noise_depth_mix_ratio) {
+  float one_minus_thickness_ratio = 1.0f - noise_thickness_mix_ratio;
+  float one_minus_depth_ratio     = 1.0f - noise_depth_mix_ratio;
+  float* tmp_thickness            = thickness_map_p;
+  float* tmp_depth                = depth_map_p;
+  float* tmp_noise                = noise_map_p;
+
+  for (int j = 0; j < dim.ly; j++) {
+    for (int i = 0; i < dim.lx;
+         i++, tmp_thickness++, tmp_depth++, tmp_noise++) {
+      *tmp_thickness = *tmp_noise * noise_thickness_mix_ratio +
+                       *tmp_thickness * one_minus_thickness_ratio;
+      *tmp_depth = *tmp_noise * noise_depth_mix_ratio +
+                   *tmp_depth * one_minus_depth_ratio;
+    }
+  }
+}
+//------------------------------------
+
+void Iwa_SoapBubbleFx::do_distance_transform(float* dst_p, TDimensionI dim,
+                                             double frame) {
+  float ar = (float)m_shape_aspect_ratio->getValue(frame);
+
+  float* f = new float[std::max(dim.lx, dim.ly)];
+
+  float max_val = 0.0f;
+
+  float* tmp_dst = dst_p;
+  /* transform along rows */
+  for (int j = 0; j < dim.ly; j++) {
+    for (int i = 0; i < dim.lx; i++, *tmp_dst++) {
+      f[i] = *tmp_dst;
+    }
+
+    tmp_dst -= dim.lx;
+
+    float* d = dt(f, dim.lx);
+    for (int i = 0; i < dim.lx; i++, tmp_dst++) {
+      *tmp_dst = d[i];
+    }
+    delete[] d;
+  }
+  /* transform along columns */
+  for (int i = 0; i < dim.lx; i++) {
+    for (int j = 0; j < dim.ly; j++) {
+      f[j] = dst_p[j * dim.lx + i];
+    }
+    float* d =
+        dt(f, dim.ly,
+           ar); /* ar : taking account of the aspect ratio of the shape */
+    for (int j = 0; j < dim.ly; j++) {
+      dst_p[j * dim.lx + i]       = d[j];
+      if (d[j] > max_val) max_val = d[j];
+    }
+    delete[] d;
+  }
+
+  tmp_dst = dst_p;
+  max_val = std::sqrt(max_val);
+
+  /* square root and normalize */
+  for (int i = 0; i < dim.lx * dim.ly; i++, *tmp_dst++) {
+    *tmp_dst = std::sqrt(*tmp_dst) / max_val;
+  }
+}
+
+//==============================================================================
+
+FX_PLUGIN_IDENTIFIER(Iwa_SoapBubbleFx, "iwa_SoapBubbleFx");
\ No newline at end of file
diff --git a/toonz/sources/stdfx/iwa_soapbubblefx.h b/toonz/sources/stdfx/iwa_soapbubblefx.h
new file mode 100644
index 0000000..2b7acb1
--- /dev/null
+++ b/toonz/sources/stdfx/iwa_soapbubblefx.h
@@ -0,0 +1,88 @@
+#pragma once
+
+/*------------------------------------
+Iwa_SoapBubbleFx
+Generates thin film interference colors from two reference images;
+one is for thickness and the other one is for shape or normal vector
+distribution of the film.
+Inherits Iwa_SpectrumFx.
+------------------------------------*/
+
+#ifndef IWA_SOAPBUBBLE_H
+#define IWA_SOAPBUBBLE_H
+
+#include "iwa_spectrumfx.h"
+
+class Iwa_SoapBubbleFx final : public Iwa_SpectrumFx {
+  FX_PLUGIN_DECLARATION(Iwa_SoapBubbleFx)
+
+protected:
+  /* target shape, used to create a pseudo normal vector */
+  TRasterFxPort m_shape;
+  /* another option, to input a depth map directly */
+  TRasterFxPort m_depth;
+  // shape parameters
+  TDoubleParamP m_binarize_threshold;
+  TDoubleParamP m_shape_aspect_ratio;
+  TDoubleParamP m_blur_radius;
+  TDoubleParamP m_blur_power;
+
+  // noise parameters
+  TIntParamP m_normal_sample_distance;
+  TIntParamP m_noise_sub_depth;
+  TDoubleParamP m_noise_resolution_s;
+  TDoubleParamP m_noise_resolution_t;
+  TDoubleParamP m_noise_sub_composite_ratio;
+  TDoubleParamP m_noise_evolution;
+  TDoubleParamP m_noise_depth_mix_ratio;
+  TDoubleParamP m_noise_thickness_mix_ratio;
+
+  template <typename RASTER, typename PIXEL>
+  void convertToBrightness(const RASTER srcRas, float* dst, TDimensionI dim);
+
+  template <typename RASTER, typename PIXEL>
+  void convertToRaster(const RASTER ras, float* thickness_map_p,
+                       float* depth_map_p, TDimensionI dim,
+                       float3* bubbleColor_p);
+
+  void processShape(double frame, TTile& shape_tile, float* depth_map_p,
+                    TDimensionI dim, const TRenderSettings& settings);
+
+  void do_binarize(TRaster32P srcRas, char* dst_p, float thres,
+                   float* distance_p, TDimensionI dim);
+
+  void do_createBlurFilter(float* dst_p, int size, float radius);
+
+  void do_applyFilter(float* depth_map_p, TDimensionI dim, float* distace_p,
+                      char* binarized_p, float* blur_filter_p,
+                      int blur_filter_size, double frame);
+
+  void processNoise(float* thickness_map_p, float* depth_map_p, TDimensionI dim,
+                    double frame, const TRenderSettings& settings);
+
+  void calc_norm_angle(float* norm_angle_p, float* depth_map_p, TDimensionI dim,
+                       int shrink);
+
+  void make_noise_map(float* noise_map_p, float* depth_map_p,
+                      float* norm_angle_p, TDimensionI dim,
+                      const QList<int>& noise_amount,
+                      const QList<QSize>& noise_base_resolution,
+                      int noise_sub_depth, float* noise_base);
+
+  float noise_interp(int left, int right, int bottom, int top, float ratio_s,
+                     float ratio_t, float* noise_layer_base, int noise_dim_x);
+
+  void add_noise(float* thickness_map_p, float* depth_map_p, TDimensionI dim,
+                 float* noise_map_p, float noise_thickness_mix_ratio,
+                 float noise_depth_mix_ratio);
+
+  void do_distance_transform(float* dst_p, TDimensionI dim, double frame);
+
+public:
+  Iwa_SoapBubbleFx();
+
+  void doCompute(TTile& tile, double frame,
+                 const TRenderSettings& settings) override;
+};
+
+#endif
\ No newline at end of file
diff --git a/toonz/sources/stdfx/iwa_spectrumfx.cpp b/toonz/sources/stdfx/iwa_spectrumfx.cpp
index 1dec6ee..d743405 100644
--- a/toonz/sources/stdfx/iwa_spectrumfx.cpp
+++ b/toonz/sources/stdfx/iwa_spectrumfx.cpp
@@ -13,25 +13,27 @@ const float PI = 3.14159265f;
 }
 
 /*------------------------------------
- シャボン色マップの生成
+ Calculate soap bubble color map
 ------------------------------------*/
-void Iwa_SpectrumFx::calcBubbleMap(float3 *bubbleColor, double frame) {
-  int j, k;     /*- bubbleColor[j][k] = [256][3] -*/
-  float d;      /*- 膜厚(μm) -*/
-  int ram;      /*- 波長のfor文用 -*/
-  float rambda; /*- 波長(μm) -*/
+void Iwa_SpectrumFx::calcBubbleMap(float3 *bubbleColor, double frame,
+                                   bool computeAngularAxis) {
+  int i, j, k;  /* bubbleColor[j][k] = [256][3] */
+  float d;      /* Thickness of the film (μm) */
+  int ram;      /* rambda iterator */
+  float rambda; /* wavelength of light (μm) */
   struct REFLECTIVITY {
-    float r_ab, t_ab, r_ba, t_ba; /*- 各境界での振幅反射率、振幅透過率 -*/
-    float r_real, r_img; /*- 薄膜の振幅反射率 -*/
-    float R;             /*- エネルギー反射率 -*/
+    /* transmission and reflection amplitudes for each boundary */
+    float r_ab, t_ab, r_ba, t_ba;
+    float r_real, r_img; /* reflection amplitude of the film */
+    float R;             /* energy reflectance */
   } p, s;
-  float R_final;                   /*- エネルギー反射率の最終版 -*/
-  float phi;                       /*- 位相 -*/
-  float color_x, color_y, color_z; /*- xyz表色系 -*/
+  float R_final;                   /* combined energy reflectance */
+  float phi;                       /* phase */
+  float color_x, color_y, color_z; /* xyz color channels */
 
   float temp_rgb_f[3];
 
-  /*- パラメータを得る -*/
+  /* obtain parameters */
   float intensity       = (float)m_intensity->getValue(frame);
   float refractiveIndex = (float)m_refractiveIndex->getValue(frame);
   float thickMax        = (float)m_thickMax->getValue(frame);
@@ -41,79 +43,99 @@ void Iwa_SpectrumFx::calcBubbleMap(float3 *bubbleColor, double frame) {
                        (float)m_BGamma->getValue(frame)};
   float lensFactor = (float)m_lensFactor->getValue(frame);
 
-  /*- 入射角は0で固定 -*/
-
-  /*- 各境界での振幅反射率、振幅透過率の計算(PS偏光とも) -*/
-  /*- P偏光 -*/
-  p.r_ab = (1.0 - refractiveIndex) / (1.0 + refractiveIndex);
-  p.t_ab = (1.0f - p.r_ab) / refractiveIndex;
-  p.r_ba = -p.r_ab;
-  p.t_ba = (1.0f + p.r_ab) * refractiveIndex;
-  /*- S偏光 -*/
-  s.r_ab = (1.0 - refractiveIndex) / (1.0 + refractiveIndex);
-  s.t_ab = 1.0f + s.r_ab;
-  s.r_ba = -s.r_ab;
-  s.t_ba = 1.0f - s.r_ab;
-
-  for (j = 0; j < 256; j++) { /*- 膜厚d -*/
-    /*- 膜厚d(μm)の計算 -*/
-    d = thickMin +
-        (thickMax - thickMin) * powf(((float)j / 255.0f), lensFactor);
-
-    /*-  膜厚が負になることもありうる。その場合は d = 0 に合わせる -*/
-    if (d < 0.0f) d = 0.0f;
-
-    /*- これから積算するので、XYZ表色系各チャンネルの初期化 -*/
-    color_x = 0.0f;
-    color_y = 0.0f;
-    color_z = 0.0f;
-
-    for (ram = 0; ram < 34; ram++) { /*- 波長λ(380nm-710nm) -*/
-      /*- 波長λ(μm)の計算 -*/
-      rambda = 0.38f + 0.01f * (float)ram;
-      /*- 位相の計算 -*/
-      phi = 4.0f * PI * refractiveIndex * d / rambda;
-      /*- 薄膜の振幅反射率の計算(PS偏光とも) -*/
-      /*- P偏光 -*/
-      p.r_real = p.r_ab + p.t_ab * p.r_ba * p.t_ba * cosf(phi);
-      p.r_img  = p.t_ab * p.r_ba * p.t_ba * sinf(phi);
-      /*- S偏光 -*/
-      s.r_real = s.r_ab + s.t_ab * s.r_ba * s.t_ba * cosf(phi);
-      s.r_img  = s.t_ab * s.r_ba * s.t_ba * sinf(phi);
-
-      p.R = p.r_real * p.r_real + p.r_img * p.r_img;
-      s.R = s.r_real * s.r_real + s.r_img * s.r_img;
-
-      /*- エネルギー反射率 -*/
-      R_final = (p.R + s.R) / 2.0f;
-
-      color_x += intensity * cie_d65[ram] * R_final * xyz[ram * 3 + 0];
-      color_y += intensity * cie_d65[ram] * R_final * xyz[ram * 3 + 1];
-      color_z += intensity * cie_d65[ram] * R_final * xyz[ram * 3 + 2];
-
-    } /*- 次のramへ(波長λ) -*/
-
-    temp_rgb_f[0] =
-        3.240479f * color_x - 1.537150f * color_y - 0.498535f * color_z;
-    temp_rgb_f[1] =
-        -0.969256f * color_x + 1.875992f * color_y + 0.041556f * color_z;
-    temp_rgb_f[2] =
-        0.055648f * color_x - 0.204043f * color_y + 1.057311f * color_z;
-
-    /*- オーバーフローをまるめる -*/
-    for (k = 0; k < 3; k++) {
-      if (temp_rgb_f[k] < 0.0f) temp_rgb_f[k] = 0.0f;
-
-      /*- ガンマ処理 -*/
-      temp_rgb_f[k] = powf((temp_rgb_f[k] / 255.0f), rgbGamma[k]);
-
-      if (temp_rgb_f[k] >= 1.0f) temp_rgb_f[k] = 1.0f;
-    }
-    bubbleColor[j].x = temp_rgb_f[0];
-    bubbleColor[j].y = temp_rgb_f[1];
-    bubbleColor[j].z = temp_rgb_f[2];
+  /* for Iwa_SpectrumFx, incident angle is fixed to 0,
+     for Iwa_SoapBubbleFx, compute for all discrete incident angles*/
+  int i_max        = (computeAngularAxis) ? 256 : 1;
+  float3 *bubble_p = bubbleColor;
+
+  /* for each discrete incident angle */
+  for (i = 0; i < i_max; i++) {
+    /* incident angle (radian) */
+    float angle_in = PI / 2.0f / 255.0f * (float)i;
+    /* refraction angle (radian) */
+    float angle_re = asinf(sinf(angle_in) / refractiveIndex);
+
+    /* transmission and reflection amplitudes for each boundary, for each
+     * polarization */
+    float cos_in = cosf(angle_in);
+    float cos_re = cosf(angle_re);
+    // P-polarized light
+    p.r_ab = (cos_re - refractiveIndex * cos_in) /
+             (cos_re + refractiveIndex * cos_re);
+    p.t_ab = (1.0f - p.r_ab) / refractiveIndex;
+    p.r_ba = -p.r_ab;
+    p.t_ba = (1.0f + p.r_ab) * refractiveIndex;
+    // S-polarized light
+    s.r_ab = (cos_in - refractiveIndex * cos_re) /
+             (cos_in + refractiveIndex * cos_re);
+    s.t_ab = 1.0f + s.r_ab;
+    s.r_ba = -s.r_ab;
+    s.t_ba = 1.0f - s.r_ab;
+
+    /* for each discrete thickness */
+    for (j = 0; j < 256; j++) {
+      /* calculate the thickness of film (μm) */
+      d = thickMin +
+          (thickMax - thickMin) * powf(((float)j / 255.0f), lensFactor);
+
+      /* there may be a case that the thickness is smaller than 0 */
+      if (d < 0.0f) d = 0.0f;
+
+      /* initialize XYZ color channels */
+      color_x = 0.0f;
+      color_y = 0.0f;
+      color_z = 0.0f;
+
+      /* for each wavelength (in the range of visible light, 380nm-710nm) */
+      for (ram = 0; ram < 34; ram++) {
+        /* wavelength `λ` (μm) */
+        rambda = 0.38f + 0.01f * (float)ram;
+        /* phase of light */
+        phi = 4.0f * PI * refractiveIndex * d * cos_re / rambda;
+        /* reflection amplitude of the film for each polarization */
+        // P-polarized light
+        p.r_real = p.r_ab + p.t_ab * p.r_ba * p.t_ba * cosf(phi);
+        p.r_img  = p.t_ab * p.r_ba * p.t_ba * sinf(phi);
+        // S-polarized light
+        s.r_real = s.r_ab + s.t_ab * s.r_ba * s.t_ba * cosf(phi);
+        s.r_img  = s.t_ab * s.r_ba * s.t_ba * sinf(phi);
+
+        p.R = p.r_real * p.r_real + p.r_img * p.r_img;
+        s.R = s.r_real * s.r_real + s.r_img * s.r_img;
+
+        /* combined energy reflectance */
+        R_final = (p.R + s.R) / 2.0f;
+
+        /* accumulate XYZ channel values */
+        color_x += intensity * cie_d65[ram] * R_final * xyz[ram * 3 + 0];
+        color_y += intensity * cie_d65[ram] * R_final * xyz[ram * 3 + 1];
+        color_z += intensity * cie_d65[ram] * R_final * xyz[ram * 3 + 2];
+
+      } /* next wavelength (ram) */
+
+      temp_rgb_f[0] =
+          3.240479f * color_x - 1.537150f * color_y - 0.498535f * color_z;
+      temp_rgb_f[1] =
+          -0.969256f * color_x + 1.875992f * color_y + 0.041556f * color_z;
+      temp_rgb_f[2] =
+          0.055648f * color_x - 0.204043f * color_y + 1.057311f * color_z;
+
+      /* clamp overflows */
+      for (k = 0; k < 3; k++) {
+        if (temp_rgb_f[k] < 0.0f) temp_rgb_f[k] = 0.0f;
+
+        /* gamma adjustment */
+        temp_rgb_f[k] = powf((temp_rgb_f[k] / 255.0f), rgbGamma[k]);
+
+        if (temp_rgb_f[k] >= 1.0f) temp_rgb_f[k] = 1.0f;
+      }
+      bubble_p->x = temp_rgb_f[0];
+      bubble_p->y = temp_rgb_f[1];
+      bubble_p->z = temp_rgb_f[2];
+      bubble_p++;
 
-  } /*- 次のjへ(膜厚d) -*/
+    } /*- next thickness d (j) -*/
+  }   /*- next incident angle (i) -*/
 }
 
 //------------------------------------
diff --git a/toonz/sources/stdfx/iwa_spectrumfx.h b/toonz/sources/stdfx/iwa_spectrumfx.h
index 6d73b75..e0ddf70 100644
--- a/toonz/sources/stdfx/iwa_spectrumfx.h
+++ b/toonz/sources/stdfx/iwa_spectrumfx.h
@@ -14,12 +14,14 @@
 
 struct float3 {
   float x, y, z;
+  float3 operator*(const float &a) { return {x * a, y * a, z * a}; }
+  float3 operator+(const float3 &a) { return {x + a.x, y + a.y, z + a.z}; }
 };
 struct float4 {
   float x, y, z, w;
 };
 
-class Iwa_SpectrumFx final : public TStandardRasterFx {
+class Iwa_SpectrumFx : public TStandardRasterFx {
   FX_PLUGIN_DECLARATION(Iwa_SpectrumFx)
 
 protected:
@@ -38,7 +40,8 @@ protected:
   TDoubleParamP m_lightIntensity;
 
   /*- シャボン色マップの生成 -*/
-  void calcBubbleMap(float3 *bubbleColor, double frame);
+  void calcBubbleMap(float3 *bubbleColor, double frame,
+                     bool computeAngularAxis = false);
 
   template <typename RASTER, typename PIXEL>
   void convertRaster(const RASTER ras, TDimensionI dim, float3 *bubbleColor);