diff --git a/apps/vaporgui/CMakeLists.txt b/apps/vaporgui/CMakeLists.txt
index 002f9b3e04..93c5bc3bed 100644
--- a/apps/vaporgui/CMakeLists.txt
+++ b/apps/vaporgui/CMakeLists.txt
@@ -237,6 +237,8 @@ set (SRCS
 	PButton.h
 	PFlowRakeRegionSelector.cpp
 	PFlowRakeRegionSelector.h
+    PFlowIntegrationRegionSelector.cpp
+    PFlowIntegrationRegionSelector.h
 	PMultiVarSelector.cpp
 	PMultiVarSelector.h
 	AppSettingsMenu.cpp
diff --git a/apps/vaporgui/FlowEventRouter.cpp b/apps/vaporgui/FlowEventRouter.cpp
index 754b4dbb66..a2379d68b1 100644
--- a/apps/vaporgui/FlowEventRouter.cpp
+++ b/apps/vaporgui/FlowEventRouter.cpp
@@ -2,6 +2,7 @@
 #include "vapor/FlowParams.h"
 #include "PWidgets.h"
 #include "PFlowRakeRegionSelector.h"
+#include "PFlowIntegrationRegionSelector.h"
 #include "PMultiVarSelector.h"
 #include "PConstantColorWidget.h"
 
@@ -10,6 +11,7 @@ typedef FlowParams FP;
 
 static RenderEventRouterRegistrar<FlowEventRouter> registrar(FlowEventRouter::GetClassType());
 const string                                       FlowEventRouter::SeedingTabName = "Seeding";
+const string                                       FlowEventRouter::IntegrationTabName = "Integration";
 
 FlowEventRouter::FlowEventRouter(QWidget *parent, ControlExec *ce) : RenderEventRouterGUI(ce, FlowParams::GetClassType())
 {
@@ -88,7 +90,7 @@ FlowEventRouter::FlowEventRouter(QWidget *parent, ControlExec *ce) : RenderEvent
     }));
     
     AddAppearanceSubtab(new PGroup({
-        (new PTFEditor(RenderParams::_colorMapVariableNameTag))->ShowOpacityBasedOnParam("NULL", 1),
+        (new PTFEditor(RenderParams::_colorMapVariableNameTag)),
         new PSection("Appearance", {
             new PConstantColorWidget,
             new PEnumDropdown(FP::RenderTypeTag, {"Tubes", "Samples", "KLGWTH"}, {FP::RenderTypeStream, FP::RenderTypeSamples, FP::RenderTypeDensity}, "Render Type"),
@@ -127,6 +129,20 @@ FlowEventRouter::FlowEventRouter(QWidget *parent, ControlExec *ce) : RenderEvent
         }))->SetTooltip("Only accessible in debug build."),
 #endif
     }));
+    
+    AddSubtab(IntegrationTabName, new PGroup({
+        new PSection("Integration Settings", {
+            new PCheckbox(FP::_doIntegrationTag, "Integrate particle values along trajectory"),
+            new PCheckbox(FP::_integrationSetAllToFinalValueTag, "Set entire stream value to integrated total"),
+            (new PVariableSelector(RenderParams::_colorMapVariableNameTag, "Scalar to Integrate"))->EnableBasedOnParam(FP::_doIntegrationTag),
+            (new PDoubleInput(FP::_integrationScalarTag, "Integration distance scale"))->EnableBasedOnParam(FP::_doIntegrationTag),
+        }),
+        (new PSection("Integration Region", {
+            new PFlowIntegrationRegionSelector1D(0),
+            new PFlowIntegrationRegionSelector1D(1),
+            new PFlowIntegrationRegionSelector1D(2),
+        }))->EnableBasedOnParam(FP::_doIntegrationTag),
+    }));
 
     AddGeometrySubtab(new PGeometrySubtab);
     AddAnnotationSubtab(new PAnnotationColorbarWidget);
diff --git a/apps/vaporgui/FlowEventRouter.h b/apps/vaporgui/FlowEventRouter.h
index b7c7df04db..819e81da46 100644
--- a/apps/vaporgui/FlowEventRouter.h
+++ b/apps/vaporgui/FlowEventRouter.h
@@ -17,6 +17,7 @@ class FlowEventRouter : public RenderEventRouterGUI {
 
 public:
     static const std::string SeedingTabName;
+    static const std::string IntegrationTabName;
 
     FlowEventRouter(QWidget *parent, VAPoR::ControlExec *ce);
     static string GetClassType() { return VAPoR::FlowRenderer::GetClassType(); }
diff --git a/apps/vaporgui/Histo.cpp b/apps/vaporgui/Histo.cpp
index 84ae09a344..57c7be16e1 100644
--- a/apps/vaporgui/Histo.cpp
+++ b/apps/vaporgui/Histo.cpp
@@ -59,6 +59,8 @@ void Histo::reset(int newNumBins)
         if (_above) delete[] _above;
         _below = nullptr;
         _above = nullptr;
+        _nBinsBelow = 0;
+        _nBinsAbove = 0;
     }
     for (int i = 0; i < _numBins; i++) _binArray[i] = 0;
     if (_below) memset(_below, 0, _nBinsBelow * sizeof(*_below));
@@ -78,6 +80,12 @@ void Histo::reset(int newNumBins, float mnData, float mxData)
     _range = _maxMapData - _minMapData;
 }
 
+void Histo::setBins(const vector<long> &bins)
+{
+    VAssert(bins.size() == _numBins);
+    for (int i = 0; i < _numBins; i++) _binArray[i] = bins[i];
+}
+
 void Histo::addToBin(float val)
 {
     // The additional checks below are because
@@ -134,12 +142,12 @@ int Histo::getMaxBinSizeBetweenIndices(const int start, const int end) const
 {
     int maxBin = 0;
 
-    if (start < 0)
+    if (start < 0 && _below)
         for (int i = max(0, start + _nBinsBelow); i < min(end + _nBinsBelow, _nBinsBelow); i++) maxBin = maxBin < _below[i] ? _below[i] : maxBin;
 
     for (int i = max(start, 0); i < min(end, _numBins); i++) maxBin = maxBin < _binArray[i] ? _binArray[i] : maxBin;
 
-    if (end >= _numBins)
+    if (end >= _numBins && _above)
         for (int i = max(start - _numBins, 0); i < min(end - _numBins, _nBinsAbove); i++) maxBin = maxBin < _above[i] ? _above[i] : maxBin;
 
     if (maxBin == 0)
@@ -221,10 +229,12 @@ int Histo::Populate(const std::string &varName, VAPoR::DataMgr *dm, VAPoR::Rende
     if (_below) {
         delete[] _below;
         _below = nullptr;
+        _nBinsBelow = 0;
     }
     if (_above) {
         delete[] _above;
         _above = nullptr;
+        _nBinsAbove = 0;
     }
 
     _getDataRange(varName, dm, rp, &_minData, &_maxData);
diff --git a/apps/vaporgui/Histo.h b/apps/vaporgui/Histo.h
index 84b1072775..774d394723 100644
--- a/apps/vaporgui/Histo.h
+++ b/apps/vaporgui/Histo.h
@@ -37,6 +37,7 @@ class Histo {
     void  reset(int newNumBins = -1);
     void  reset(int newNumBins, float mnData, float mxData);
     void  addToBin(float val);
+    void  setBins(const vector<long> &bins);
     int   getMaxBinSize();
     int   getMaxBinSizeBetweenIndices(const int start, const int end) const;
     int   getNumBins() const;
diff --git a/apps/vaporgui/PFlowIntegrationRegionSelector.cpp b/apps/vaporgui/PFlowIntegrationRegionSelector.cpp
new file mode 100644
index 0000000000..f4e7da9033
--- /dev/null
+++ b/apps/vaporgui/PFlowIntegrationRegionSelector.cpp
@@ -0,0 +1,6 @@
+#include "PFlowIntegrationRegionSelector.h"
+#include <vapor/FlowParams.h>
+
+using VAPoR::FlowParams;
+
+VAPoR::Box *PFlowIntegrationRegionSelector1D::getBox() const { return getParams<FlowParams>()->GetIntegrationBox(); }
diff --git a/apps/vaporgui/PFlowIntegrationRegionSelector.h b/apps/vaporgui/PFlowIntegrationRegionSelector.h
new file mode 100644
index 0000000000..85d32e32c8
--- /dev/null
+++ b/apps/vaporgui/PFlowIntegrationRegionSelector.h
@@ -0,0 +1,11 @@
+#pragma once
+
+#include "PRegionSelector.h"
+
+class PFlowIntegrationRegionSelector1D : public PRegionSelector1D {
+public:
+    using PRegionSelector1D::PRegionSelector1D;
+
+protected:
+    VAPoR::Box *getBox() const override;
+};
diff --git a/apps/vaporgui/TFHistogramWidget.cpp b/apps/vaporgui/TFHistogramWidget.cpp
index c7525b15fa..1d024920df 100644
--- a/apps/vaporgui/TFHistogramWidget.cpp
+++ b/apps/vaporgui/TFHistogramWidget.cpp
@@ -42,9 +42,18 @@ void TFHistogramMap::PopulateSettingsMenu(QMenu *menu) const
 
 void TFHistogramMap::paramsUpdate()
 {
-    int numBins = width() ? width() : 256;
-    if (_histo.getNumBins() != numBins) _histo.reset(numBins);    // This can be called before it is resized
-    if (_histo.PopulateIfNeeded(getVariableName(), getDataMgr(), getRenderParams()) < 0) MSG_ERR("Failed to populate histogram");
+    RenderParams *rp = getRenderParams();
+    if (rp->GetValueDoubleVec(RenderParams::CustomHistogramRangeTag).size()) {
+        auto range = rp->GetValueDoubleVec(RenderParams::CustomHistogramRangeTag);
+        auto bins = rp->GetValueLongVec(RenderParams::CustomHistogramDataTag);
+        int  numBins = bins.size();
+        _histo.reset(numBins, range[0], range[1]);
+        _histo.setBins(bins);
+    } else {
+        int numBins = width() ? width() : 256;
+        if (_histo.getNumBins() != numBins) _histo.reset(numBins);    // This can be called before it is resized
+        if (_histo.PopulateIfNeeded(getVariableName(), getDataMgr(), getRenderParams()) < 0) MSG_ERR("Failed to populate histogram");
+    }
 
     _scalingMenu->Update(getRenderParams());
     update();
diff --git a/apps/vaporgui/TFHistogramWidget.h b/apps/vaporgui/TFHistogramWidget.h
index 9ac41bf692..377de22850 100644
--- a/apps/vaporgui/TFHistogramWidget.h
+++ b/apps/vaporgui/TFHistogramWidget.h
@@ -30,8 +30,6 @@ class TFHistogramMap : public TFMap {
     //    void mouseDoubleClickEvent(QMouseEvent *event);
 
 private:
-    VAPoR::DataMgr *        _dataMgr = nullptr;
-    VAPoR::RenderParams *   _renderParams = nullptr;
     Histo                   _histo;
     ParamsDropdownMenuItem *_scalingMenu;
     bool                    _dynamicScaling = true;
diff --git a/apps/vaporgui/TFMappingRangeSelector.cpp b/apps/vaporgui/TFMappingRangeSelector.cpp
index b5ea89a994..6c066c0769 100644
--- a/apps/vaporgui/TFMappingRangeSelector.cpp
+++ b/apps/vaporgui/TFMappingRangeSelector.cpp
@@ -3,6 +3,8 @@
 #include <vapor/ParamsMgr.h>
 #include <vapor/DataMgrUtils.h>
 
+using VAPoR::RenderParams;
+
 TFMappingRangeSelector::TFMappingRangeSelector(const std::string &variableNameTag) : _variableNameTag(variableNameTag)
 {
     AllowCustomRange();
@@ -27,7 +29,7 @@ void TFMappingRangeSelector::Update(VAPoR::DataMgr *dataMgr, VAPoR::ParamsMgr *p
         min = range[0];
         max = range[1];
     } else {
-        _getDataRange(dataMgr, rParams, &min, &max);
+        _getDefaultRange(dataMgr, rParams, &min, &max);
     }
     SetRange(min, max);
 
@@ -35,6 +37,17 @@ void TFMappingRangeSelector::Update(VAPoR::DataMgr *dataMgr, VAPoR::ParamsMgr *p
     SetValue(mapperRange[0], mapperRange[1]);
 }
 
+void TFMappingRangeSelector::_getDefaultRange(VAPoR::DataMgr *d, VAPoR::RenderParams *r, float *min, float *max) const
+{
+    if (r->GetValueDoubleVec(RenderParams::CustomHistogramRangeTag).size() == 2) {
+        auto range = r->GetValueDoubleVec(RenderParams::CustomHistogramRangeTag);
+        *min = range[0];
+        *max = range[1];
+    } else {
+        _getDataRange(d, r, min, max);
+    }
+}
+
 void TFMappingRangeSelector::_getDataRange(VAPoR::DataMgr *d, VAPoR::RenderParams *r, float *min, float *max) const
 {
     vector<double> minExt, maxExt;
diff --git a/apps/vaporgui/TFMappingRangeSelector.h b/apps/vaporgui/TFMappingRangeSelector.h
index 3e5d32fc93..e1ca40318c 100644
--- a/apps/vaporgui/TFMappingRangeSelector.h
+++ b/apps/vaporgui/TFMappingRangeSelector.h
@@ -22,6 +22,7 @@ class TFMappingRangeSelector : public QRangeSliderTextCombo {
     VAPoR::ParamsMgr *   _paramsMgr = nullptr;
     const std::string &  _variableNameTag;
 
+    void                   _getDefaultRange(VAPoR::DataMgr *dataMgr, VAPoR::RenderParams *rParams, float *min, float *max) const;
     void                   _getDataRange(VAPoR::DataMgr *dataMgr, VAPoR::RenderParams *rParams, float *min, float *max) const;
     std::string            _getVariableName() const;
     VAPoR::MapperFunction *_getTF() const;
diff --git a/apps/vaporgui/VizWin.cpp b/apps/vaporgui/VizWin.cpp
index 4a5ebe24ac..058f05b307 100644
--- a/apps/vaporgui/VizWin.cpp
+++ b/apps/vaporgui/VizWin.cpp
@@ -545,7 +545,7 @@ string VizWin::_getCurrentMouseMode() const
     GUIStateParams *guiP = (GUIStateParams *)paramsMgr->GetParams(GUIStateParams::GetClassType());
 
     string activeTab = guiP->ActiveTab();
-    if (activeTab == RenderEventRouterGUI::GeometryTabName || activeTab == FlowEventRouter::SeedingTabName)
+    if (activeTab == RenderEventRouterGUI::GeometryTabName || activeTab == FlowEventRouter::SeedingTabName || activeTab == FlowEventRouter::IntegrationTabName)
         return MouseModeParams::GetRegionModeName();
     else
         return MouseModeParams::GetNavigateModeName();
@@ -585,6 +585,10 @@ void VizWin::_setNewExtents()
                 fp->SetRake(b);
             }
         }
+    } else if (_manipFlowIntegrationFlag) {
+        FlowParams *fp = dynamic_cast<FlowParams *>(rParams);
+        VAssert(fp);
+        fp->GetIntegrationBox()->SetExtents(llc, urc);
     } else {
         box->SetExtents(llc, urc);
     }
@@ -814,10 +818,12 @@ void VizWin::updateManip(bool initialize)
         urc = maxExts;
     } else {
         _manipFlowSeedFlag = false;
+        _manipFlowIntegrationFlag = false;
         GUIStateParams *gp = (GUIStateParams *)_controlExec->GetParamsMgr()->GetParams(GUIStateParams::GetClassType());
 
         if (rParams->GetName() == FlowParams::GetClassType()) {
             if (gp->ActiveTab() == FlowEventRouter::SeedingTabName) _manipFlowSeedFlag = true;
+            if (gp->ActiveTab() == FlowEventRouter::IntegrationTabName) _manipFlowIntegrationFlag = true;
         }
         if (_manipFlowSeedFlag) {
             FlowParams *fp = dynamic_cast<FlowParams *>(rParams);
@@ -853,6 +859,10 @@ void VizWin::updateManip(bool initialize)
             urc[1] = b[3];
             llc[2] = b[4];
             urc[2] = b[5];
+        } else if (_manipFlowIntegrationFlag) {
+            FlowParams *fp = dynamic_cast<FlowParams *>(rParams);
+            VAssert(fp);
+            fp->GetIntegrationBox()->GetExtents(llc, urc);
         } else {
             VAPoR::Box *box = rParams->GetBox();
             box->GetExtents(llc, urc);
diff --git a/apps/vaporgui/VizWin.h b/apps/vaporgui/VizWin.h
index b74f824a47..bcc0fcc73a 100644
--- a/apps/vaporgui/VizWin.h
+++ b/apps/vaporgui/VizWin.h
@@ -132,6 +132,7 @@ public slots:
     bool       _navigateFlag;
     bool       _manipFlag;
     bool       _manipFlowSeedFlag = false;
+    bool       _manipFlowIntegrationFlag = false;
     Trackball *_trackBall;
 
     std::vector<double> _getScreenCoords(QMouseEvent *e) const;
diff --git a/include/vapor/Advection.h b/include/vapor/Advection.h
index 10c94eb6da..5756422544 100644
--- a/include/vapor/Advection.h
+++ b/include/vapor/Advection.h
@@ -38,8 +38,13 @@ class FLOW_API Advection final {
     //   If "skipNonZero" is true, then this function only overwrites zeros.
     //   Otherwise, it will overwrite values anyway.
     int CalculateParticleValues(Field *scalarField, bool skipNonZero);
+    int CalculateParticleIntegratedValues(Field *scalarField, const bool skipNonZero, const float distScale = 1.f, const std::vector<double> &integrateWithinVolumeMin = {-FLT_MAX, -FLT_MAX, -FLT_MAX},
+                                          const std::vector<double> &integrateWithinVolumeMax = {FLT_MAX, FLT_MAX, FLT_MAX});
+    void SetAllStreamValuesToFinalValue(int realNSamples);
     int CalculateParticleProperties(Field *scalarField);
 
+    void CalculateParticleHistogram(std::vector<double> &bounds, std::vector<long> &bins);
+
     // Reset all particle values to zero
     void ResetParticleValues();
     // Clear all existing properties of a particle
@@ -100,6 +105,10 @@ class FLOW_API Advection final {
     // Adjust input "val" according to the bound specified by min and max.
     // Returns the value after adjustment.
     float _applyPeriodic(float val, float min, float max) const;
+
+    void        _calculateParticleIntegratedValue(Particle &p, const Particle &prev, const Field *scalarField, const bool skipNonZero, const float distScale,
+                                                  const std::vector<double> &integrateWithinVolumeMin, const std::vector<double> &integrateWithinVolumeMax) const;
+    static bool _isParticleInsideVolume(const Particle &p, const std::vector<double> &min, const std::vector<double> &max);
 };
 };    // namespace flow
 
diff --git a/include/vapor/FlowParams.h b/include/vapor/FlowParams.h
index 5ae3d26c43..ec94a72d42 100644
--- a/include/vapor/FlowParams.h
+++ b/include/vapor/FlowParams.h
@@ -17,15 +17,20 @@ enum class FlowDir : int { FORWARD = 0, BACKWARD = 1, BI_DIR = 2 };
 
 class FlowParams;
 class PARAMS_API FakeRakeBox : public Box {
+    string _tag;
+
 public:
     using Box::Box;
     FlowParams *parent = nullptr;
+
+    void        Initialize(string tag);
     void        SetExtents(const vector<double> &minExt, const vector<double> &maxExt) override;
 };
 
 class PARAMS_API FlowParams : public RenderParams {
     StateSave    _fakeRakeStateSave;
     FakeRakeBox *_fakeRakeBox = nullptr;
+    FakeRakeBox *_fakeIntegrationBox = nullptr;
     bool         _initialized = false;
 
 public:
@@ -86,6 +91,9 @@ class PARAMS_API FlowParams : public RenderParams {
     std::vector<float> GetRake() const;
     void               SetRake(const std::vector<float> &);
 
+    Box *GetIntegrationBox();
+    void SetIntegrationVolume(const std::vector<float> &);
+
     /*
      *This result vector could be of size 2 or 3.
      */
@@ -170,6 +178,10 @@ class PARAMS_API FlowParams : public RenderParams {
     static const std::string _yPeriodicTag;
     static const std::string _zPeriodicTag;
     static const std::string _rakeTag;
+    static const std::string _doIntegrationTag;
+    static const std::string _integrationScalarTag;
+    static const std::string _integrationSetAllToFinalValueTag;
+    static const std::string _integrationBoxTag;
     static const std::string _rakeBiasVariable;
     static const std::string _rakeBiasStrength;
     static const std::string _pastNumOfTimeSteps;
diff --git a/include/vapor/FlowRenderer.h b/include/vapor/FlowRenderer.h
index b9539a01c8..9eb2c24815 100644
--- a/include/vapor/FlowRenderer.h
+++ b/include/vapor/FlowRenderer.h
@@ -76,6 +76,10 @@ class RENDER_API FlowRenderer final : public Renderer {
     FlowStatus         _renderStatus = FlowStatus::SIMPLE_OUTOFDATE;
     std::string        _cache_rakeBiasVariable;
     std::string        _cache_seedInputFilename;
+    bool                _cache_doIntegration;
+    bool                _cache_integrationSetAllToFinalValue;
+    float               _cache_integrationDistScalar;
+    std::vector<double> _cache_integrationVolume;
 
     // This Advection class is only used in bi-directional advection mode
     std::unique_ptr<flow::Advection> _2ndAdvection;
diff --git a/include/vapor/RenderParams.h b/include/vapor/RenderParams.h
index e4ab5e8417..cff3c5e8e2 100644
--- a/include/vapor/RenderParams.h
+++ b/include/vapor/RenderParams.h
@@ -436,6 +436,8 @@ class PARAMS_API RenderParams : public ParamsBase {
     static const string _yFieldVariableNameTag;
     static const string _zFieldVariableNameTag;
     static const string _constantOpacityTag;
+    static const string CustomHistogramDataTag;
+    static const string CustomHistogramRangeTag;
 
     //! If a renderer supports rotation about a point of origin,
     //! this string identifies the parameter for the origin's
diff --git a/lib/flow/Advection.cpp b/lib/flow/Advection.cpp
index 6dafcbb6dd..5e6469ad1b 100644
--- a/lib/flow/Advection.cpp
+++ b/lib/flow/Advection.cpp
@@ -325,6 +325,108 @@ int Advection::CalculateParticleValues(Field *scalar, bool skipNonZero)
     return 0;
 }
 
+
+int Advection::CalculateParticleIntegratedValues(Field *scalar, const bool skipNonZero, const float distScale, const std::vector<double> &integrateWithinVolumeMin,
+                                                 const std::vector<double> &integrateWithinVolumeMax)
+{
+    // For steady fields, we calculate values one stream at a time
+    if (scalar->IsSteady) {
+        if (scalar->LockParams() != 0) return PARAMS_ERROR;
+
+        _valueVarName = scalar->ScalarName;
+
+        for (auto &s : _streams) {
+            if (s.size() && !s[0].IsSpecial()) s[0].value = 0;
+
+            for (int i = 1; i < s.size(); i++) {
+                auto &prev = s[i - 1];
+                auto &p = s[i];
+                _calculateParticleIntegratedValue(p, prev, scalar, skipNonZero, distScale, integrateWithinVolumeMin, integrateWithinVolumeMax);
+            }
+        }
+
+        scalar->UnlockParams();
+    }
+    // For unsteady fields, we calculate values at one timestep at a time
+    else {
+        size_t mostSteps = 0;
+        for (auto &s : _streams) {
+            if (s.size() > mostSteps) mostSteps = s.size();
+        }
+
+        _valueVarName = scalar->ScalarName;
+
+        for (auto &s : _streams)
+            if (s.size() && !s[0].IsSpecial()) s[0].value = 0;
+
+        for (size_t i = 1; i < mostSteps; i++) {
+            for (auto &s : _streams) {
+                if (i < s.size()) {
+                    auto &prev = s[i - 1];
+                    auto &p = s[i];
+                    _calculateParticleIntegratedValue(p, prev, scalar, skipNonZero, distScale, integrateWithinVolumeMin, integrateWithinVolumeMax);
+                }
+            }    // end of a stream
+        }        // end of all steps
+    }
+
+    return 0;
+}
+
+void Advection::_calculateParticleIntegratedValue(Particle &p, const Particle &prev, const Field *scalarField, const bool skipNonZero, const float distScale,
+                                                  const std::vector<double> &integrateWithinVolumeMin, const std::vector<double> &integrateWithinVolumeMax) const
+{
+    // Skip this particle if it is a separator
+    if (p.IsSpecial()) return;
+    if (prev.IsSpecial()) {
+        p.value = 0;
+        return;
+    }
+
+    // Do not evaluate this particle if its value is non-zero
+    if (skipNonZero && p.value != 0.0f) return;
+
+    if (!_isParticleInsideVolume(p, integrateWithinVolumeMin, integrateWithinVolumeMax)) {
+        p.value = prev.value;
+        return;
+    }
+
+    float value;
+    int   rv = scalarField->GetScalar(p.time, p.location, value);
+    if (rv != 0) {    // If non-0, then outside the volume
+        p.value = prev.value;
+        return;
+    }
+
+    float dist = glm::distance(prev.location, p.location);
+    p.value = prev.value + value * dist * distScale;
+}
+
+void Advection::SetAllStreamValuesToFinalValue(int realNSamples)
+{
+    for (auto &s : _streams) {
+        float finalValue = 0;
+
+        int sampleCount = 0;
+        for (auto &p : s) {
+            if (!p.IsSpecial()) {
+                finalValue = p.value;
+                sampleCount++;
+            }
+            if (sampleCount == realNSamples) break;
+        }
+
+        int setCount = 0;
+        for (auto &p : s) {
+            if (!p.IsSpecial()) {
+                p.value = finalValue;
+                setCount++;
+            }
+            if (setCount == sampleCount) break;
+        }
+    }
+}
+
 int Advection::CalculateParticleProperties(Field *scalar)
 {
     // Test if this scalar property is already calculated.
@@ -380,6 +482,30 @@ int Advection::CalculateParticleProperties(Field *scalar)
     return 0;
 }
 
+void Advection::CalculateParticleHistogram(std::vector<double> &outBounds, std::vector<long> &bins)
+{
+    std::vector<float> samples;
+
+    for (const auto &s : _streams)
+        for (const auto &p : s)
+            if (!p.IsSpecial()) samples.push_back(p.value);
+
+    auto  bounds = std::minmax_element(samples.begin(), samples.end());
+    float minValue = *bounds.first;
+    float maxValue = *bounds.second;
+    float range = maxValue - minValue;
+
+    int nBins = bins.size();
+    assert(nBins != 0);
+    std::fill(bins.begin(), bins.end(), 0);
+
+    for (const auto &s : samples) { bins[std::min(nBins - 1, std::max(0, (int)((nBins - 1) * (s - minValue) / range)))]++; }
+
+    outBounds.resize(2);
+    outBounds[0] = minValue;
+    outBounds[1] = maxValue;
+}
+
 int Advection::_advectEuler(Field *velocity, const Particle &p0, double dt, Particle &p1) const
 {
     glm::vec3 v0;
@@ -530,3 +656,10 @@ float Advection::_applyPeriodic(float val, float min, float max) const
 auto Advection::GetValueVarName() const -> std::string { return _valueVarName; }
 
 auto Advection::GetPropertyVarNames() const -> std::vector<std::string> { return _propertyVarNames; }
+
+bool Advection::_isParticleInsideVolume(const Particle &p, const std::vector<double> &min, const std::vector<double> &max)
+{
+    if (p.location[0] < min[0] || p.location[1] < min[1] || p.location[0] > max[0] || p.location[1] > max[1]) { return false; }
+    if (min.size() > 2 && (p.location[2] < min[2] || p.location[2] > max[2])) { return false; }
+    return true;
+}
diff --git a/lib/params/FlowParams.cpp b/lib/params/FlowParams.cpp
index 61e44133ff..e144fd757b 100644
--- a/lib/params/FlowParams.cpp
+++ b/lib/params/FlowParams.cpp
@@ -37,6 +37,10 @@ const std::string FlowParams::_xPeriodicTag = "PeriodicTag_X";
 const std::string FlowParams::_yPeriodicTag = "PeriodicTag_Y";
 const std::string FlowParams::_zPeriodicTag = "PeriodicTag_Z";
 const std::string FlowParams::_rakeTag = "RakeTag";
+const std::string FlowParams::_doIntegrationTag = "DoIntegrationTag";
+const std::string FlowParams::_integrationScalarTag = "IntegrationScalarTag";
+const std::string FlowParams::_integrationSetAllToFinalValueTag = "IntegrationSetAllToFinalValueTag";
+const std::string FlowParams::_integrationBoxTag = "IntegrationBoxTag";
 const std::string FlowParams::_rakeBiasVariable = "RakeBiasVariable";
 const std::string FlowParams::_rakeBiasStrength = "RakeBiasStrength";
 const std::string FlowParams::_pastNumOfTimeSteps = "PastNumOfTimeSteps";
@@ -94,6 +98,7 @@ FlowParams::~FlowParams()
 {
     SetDiagMsg("FlowParams::~FlowParams() this=%p", this);
     if (_fakeRakeBox) delete _fakeRakeBox;
+    if (_fakeIntegrationBox) delete _fakeIntegrationBox;
 }
 
 int FlowParams::Initialize()
@@ -116,6 +121,7 @@ int FlowParams::Initialize()
         floats[i * 2 + 1] = maxext[i];
     }
     this->SetRake(floats);
+    this->SetIntegrationVolume(floats);
 
     vector<float> rake;
     rake.push_back(minext[0]);
@@ -127,6 +133,7 @@ int FlowParams::Initialize()
         rake.push_back(maxext[2]);
     }
     SetRake(rake);
+    SetIntegrationVolume(rake);
 
     SetFlowDirection((int)FlowDir::FORWARD);
     SetSteadyNumOfSteps(100);
@@ -138,6 +145,9 @@ int FlowParams::Initialize()
     SetSeedInputFilename(Wasp::GetSharePath("examples/listOfSeeds.txt"));
     SetIsSteady(true);
     SetPastNumOfTimeSteps(std::max(_dataMgr->GetNumTimeSteps() - 1, 1));
+    SetValueLong(_doIntegrationTag, "", false);
+    SetValueDouble(_integrationScalarTag, "", 1.f);
+    SetValueLong(_integrationSetAllToFinalValueTag, "", false);
 
     return (0);
 }
@@ -205,11 +215,36 @@ void FlowParams::SetPeriodic(const std::vector<bool> &bools)
     if (bools.size() == 3) SetValueLong(_zPeriodicTag, "", bools[2]);
 }
 
+void FakeRakeBox::Initialize(string tag)
+{
+    _tag = tag;
+
+    vector<double> flowBox = parent->GetValueDoubleVec(tag);
+    if (flowBox.size() != 4 && flowBox.size() != 6) {
+        vector<double> min, max;
+        parent->GetBox()->GetExtents(min, max);
+        SetExtents(min, max);
+        flowBox = parent->GetValueDoubleVec(tag);
+    }
+    vector<double> min, max;
+
+    min.push_back(flowBox[0]);
+    max.push_back(flowBox[1]);
+    min.push_back(flowBox[2]);
+    max.push_back(flowBox[3]);
+    if (flowBox.size() == 6) {
+        min.push_back(flowBox[4]);
+        max.push_back(flowBox[5]);
+    }
+
+    Box::SetExtents(min, max);
+}
+
 void FakeRakeBox::SetExtents(const vector<double> &min, const vector<double> &max)
 {
     VAssert(parent);
 
-    vector<float> rake;
+    vector<double> rake;
     rake.push_back(min[0]);
     rake.push_back(max[0]);
     rake.push_back(min[1]);
@@ -219,7 +254,7 @@ void FakeRakeBox::SetExtents(const vector<double> &min, const vector<double> &ma
         rake.push_back(max[2]);
     }
 
-    parent->SetRake(rake);
+    parent->SetValueDoubleVec(_tag, "", rake);
 }
 
 Box *FlowParams::GetRakeBox()
@@ -233,23 +268,25 @@ Box *FlowParams::GetRakeBox()
     _fakeRakeBox->SetPlanar(extBox->IsPlanar());
     _fakeRakeBox->SetOrientation(extBox->GetOrientation());
 
-    const auto rake = GetRake();
-    const auto rakesize = rake.size();
-    VAssert(rakesize == 4 || rakesize == 6);
-    vector<double> min, max;
+    _fakeRakeBox->Initialize(_rakeTag);
+
+    return _fakeRakeBox;
+}
 
-    min.push_back(rake[0]);
-    max.push_back(rake[1]);
-    min.push_back(rake[2]);
-    max.push_back(rake[3]);
-    if (rakesize == 6) {
-        min.push_back(rake[4]);
-        max.push_back(rake[5]);
+Box *FlowParams::GetIntegrationBox()
+{
+    if (!_fakeIntegrationBox) {
+        _fakeIntegrationBox = new FakeRakeBox(&_fakeRakeStateSave);
+        _fakeIntegrationBox->parent = this;
     }
 
-    _fakeRakeBox->Box::SetExtents(min, max);
+    Box *extBox = GetBox();
+    _fakeIntegrationBox->SetPlanar(extBox->IsPlanar());
+    _fakeIntegrationBox->SetOrientation(extBox->GetOrientation());
 
-    return _fakeRakeBox;
+    _fakeIntegrationBox->Initialize(_integrationBoxTag);
+
+    return _fakeIntegrationBox;
 }
 
 std::vector<float> FlowParams::GetRake() const
@@ -276,6 +313,15 @@ void FlowParams::SetRake(const std::vector<float> &rake)
     SetValueDoubleVec(_rakeTag, "rake boundaries", doubles);
 }
 
+void FlowParams::SetIntegrationVolume(const std::vector<float> &rake)
+{
+    const auto rakesize = rake.size();
+    VAssert(rakesize == 4 || rakesize == 6);
+    std::vector<double> doubles(rakesize);
+    std::copy(rake.cbegin(), rake.cend(), doubles.begin());
+    SetValueDoubleVec(_integrationBoxTag, "", doubles);
+}
+
 std::vector<long> FlowParams::GetGridNumOfSeeds() const
 {
     vector<long> sav(3);
diff --git a/lib/params/RenderParams.cpp b/lib/params/RenderParams.cpp
index 42b0d543bb..3e68826123 100644
--- a/lib/params/RenderParams.cpp
+++ b/lib/params/RenderParams.cpp
@@ -46,6 +46,8 @@ const string RenderParams::_variableNameTag = "VariableName";
 const string RenderParams::_useSingleColorTag = "UseSingleColor";
 const string RenderParams::_constantColorTag = "ConstantColor";
 const string RenderParams::_constantOpacityTag = "ConstantOpacity";
+const string RenderParams::CustomHistogramDataTag = "CustomHistogramData";
+const string RenderParams::CustomHistogramRangeTag = "CustomHistogramRange";
 const string RenderParams::_CompressionLevelTag = "CompressionLevel";
 const string RenderParams::_RefinementLevelTag = "RefinementLevel";
 const string RenderParams::_transferFunctionsTag = "MapperFunctions";
diff --git a/lib/render/FlowRenderer.cpp b/lib/render/FlowRenderer.cpp
index b13316efaf..3a5ef9381f 100644
--- a/lib/render/FlowRenderer.cpp
+++ b/lib/render/FlowRenderer.cpp
@@ -303,10 +303,45 @@ int FlowRenderer::_paintGL(bool fast)
     }
 
     if (!_coloringComplete) {
-        rv = _advection.CalculateParticleValues(&_colorField, true);
-        _printNonZero(rv, __FILE__, __func__, __LINE__);
-        if (_2ndAdvection)    // bi-directional advection
-            rv = _2ndAdvection->CalculateParticleValues(&_colorField, true);
+        bool integrate = params->GetValueLong(params->_doIntegrationTag, false);
+        bool setAllToFinalValue = params->GetValueLong(params->_integrationSetAllToFinalValueTag, false);
+
+        if (integrate) {
+            vector<double> integrationVolumeMin, integrationVolumeMax;
+            params->GetIntegrationBox()->GetExtents(integrationVolumeMin, integrationVolumeMax);
+            float distScale = params->GetValueDouble(params->_integrationScalarTag, 1.f);
+            _advection.CalculateParticleIntegratedValues(&_colorField, true, distScale, integrationVolumeMin, integrationVolumeMax);
+            _printNonZero(rv, __FILE__, __func__, __LINE__);
+            if (_2ndAdvection)    // bi-directional advection
+                rv = _2ndAdvection->CalculateParticleIntegratedValues(&_colorField, true, distScale, integrationVolumeMin, integrationVolumeMax);
+
+            if (setAllToFinalValue) {
+                int numSamplesPerStream;
+                if (_cache_isSteady)
+                    numSamplesPerStream = _cache_steadyNumOfSteps;
+                else
+                    numSamplesPerStream = _cache_currentTS;
+                _advection.SetAllStreamValuesToFinalValue(numSamplesPerStream);
+                if (_2ndAdvection) _2ndAdvection->SetAllStreamValuesToFinalValue(numSamplesPerStream);
+            }
+
+            vector<double> histoRange;
+            vector<long>   histo(256);
+            _advection.CalculateParticleHistogram(histoRange, histo);
+
+            params->SetValueLongVec(RenderParams::CustomHistogramDataTag, "", histo);
+            params->SetValueDoubleVec(RenderParams::CustomHistogramRangeTag, "", histoRange);
+        } else {
+            rv = _advection.CalculateParticleValues(&_colorField, true);
+            _printNonZero(rv, __FILE__, __func__, __LINE__);
+            if (_2ndAdvection)    // bi-directional advection
+                rv = _2ndAdvection->CalculateParticleValues(&_colorField, true);
+
+            if (params->GetValueDoubleVec(RenderParams::CustomHistogramRangeTag).size()) {
+                params->SetValueLongVec(RenderParams::CustomHistogramDataTag, "", {});
+                params->SetValueDoubleVec(RenderParams::CustomHistogramRangeTag, "", {});
+            }
+        }
         _printNonZero(rv, __FILE__, __func__, __LINE__);
         _coloringComplete = true;
     }
@@ -856,6 +891,19 @@ int FlowRenderer::_updateFlowCacheAndStates(const FlowParams *params)
         }
     }
 
+    const auto doIntegration = params->GetValueLong(params->_doIntegrationTag, false);
+    const auto integrationSetAllToFinalValue = params->GetValueLong(params->_integrationSetAllToFinalValueTag, false);
+    const auto integrationDistScalar = params->GetValueDouble(params->_integrationScalarTag, false);
+    const auto integrationVolume = params->GetValueDoubleVec(params->_integrationBoxTag);
+    if (doIntegration != _cache_doIntegration || integrationSetAllToFinalValue != _cache_integrationSetAllToFinalValue || integrationDistScalar != _cache_integrationDistScalar
+        || integrationVolume != _cache_integrationVolume) {
+        _colorStatus = FlowStatus::SIMPLE_OUTOFDATE;
+    }
+    _cache_doIntegration = doIntegration;
+    _cache_integrationSetAllToFinalValue = integrationSetAllToFinalValue;
+    _cache_integrationDistScalar = integrationDistScalar;
+    _cache_integrationVolume = integrationVolume;
+
     //
     // Now we branch into steady and unsteady cases, and treat them separately
     //
@@ -867,6 +915,7 @@ int FlowRenderer::_updateFlowCacheAndStates(const FlowParams *params)
                 if (_velocityStatus == FlowStatus::UPTODATE) _velocityStatus = FlowStatus::TIME_STEP_OOD;
             }
             if (params->GetSteadyNumOfSteps() != _cache_steadyNumOfSteps) _renderStatus = FlowStatus::SIMPLE_OUTOFDATE;
+            if (params->GetSteadyNumOfSteps() < _cache_steadyNumOfSteps && doIntegration) _colorStatus = FlowStatus::SIMPLE_OUTOFDATE;
             _cache_steadyNumOfSteps = params->GetSteadyNumOfSteps();
 
             if (_cache_currentTS != params->GetCurrentTimestep()) {
diff --git a/share/shaders/FlowGlyphsArrow.frag b/share/shaders/FlowGlyphsArrow.frag
index 518d695950..06b70820fa 100644
--- a/share/shaders/FlowGlyphsArrow.frag
+++ b/share/shaders/FlowGlyphsArrow.frag
@@ -56,4 +56,6 @@ void main() {
         color.rgb *= PhongLighting(normal, ld);
     }
     fragment = color;
+    
+    if (color.a < 0.05) discard;
 }
diff --git a/share/shaders/FlowGlyphsArrow2D.frag b/share/shaders/FlowGlyphsArrow2D.frag
index 3519b174d4..ec5bec71aa 100644
--- a/share/shaders/FlowGlyphsArrow2D.frag
+++ b/share/shaders/FlowGlyphsArrow2D.frag
@@ -36,5 +36,7 @@ void main()
 
 
         fragment = c;
+        
+        if (c.a < 0.05) discard;
     }
 }
diff --git a/share/shaders/FlowGlyphsLineDirArrow2D.frag b/share/shaders/FlowGlyphsLineDirArrow2D.frag
index 94bd8d5fe8..b260fdab4f 100644
--- a/share/shaders/FlowGlyphsLineDirArrow2D.frag
+++ b/share/shaders/FlowGlyphsLineDirArrow2D.frag
@@ -45,4 +45,6 @@ void main()
         discard;
 
     fragment = color;
+    
+    if (color.a < 0.05) discard;
 }
diff --git a/share/shaders/FlowGlyphsSphere2D.frag b/share/shaders/FlowGlyphsSphere2D.frag
index 3519b174d4..ec5bec71aa 100644
--- a/share/shaders/FlowGlyphsSphere2D.frag
+++ b/share/shaders/FlowGlyphsSphere2D.frag
@@ -36,5 +36,7 @@ void main()
 
 
         fragment = c;
+        
+        if (c.a < 0.05) discard;
     }
 }
diff --git a/share/shaders/FlowGlyphsSphereSplat.frag b/share/shaders/FlowGlyphsSphereSplat.frag
index 7cdfaab2ac..ca74e78a3c 100644
--- a/share/shaders/FlowGlyphsSphereSplat.frag
+++ b/share/shaders/FlowGlyphsSphereSplat.frag
@@ -50,4 +50,6 @@ void main() {
 
 	color.rgb *= PhongLighting(t2);
     fragment = color;
+    
+    if (color.a < 0.05) discard;
 }
diff --git a/share/shaders/FlowGlyphsTubeDirArrow.frag b/share/shaders/FlowGlyphsTubeDirArrow.frag
index 92712c49e1..975504d3c8 100644
--- a/share/shaders/FlowGlyphsTubeDirArrow.frag
+++ b/share/shaders/FlowGlyphsTubeDirArrow.frag
@@ -61,4 +61,6 @@ void main() {
         discard;
 
     fragment = color;
+    
+    if (color.a < 0.05) discard;
 }
diff --git a/share/shaders/FlowLine.frag b/share/shaders/FlowLine.frag
index e601642cfb..8d0cc4f407 100644
--- a/share/shaders/FlowLine.frag
+++ b/share/shaders/FlowLine.frag
@@ -20,6 +20,8 @@ void main(void)
     {
         float valTranslate = (scalarV - colorMapRange.x) / colorMapRange.z;
         color = texture( colorMapTexture, valTranslate );
+        
+        if (color.a < 0.05) discard;
     }
 }
 
diff --git a/share/shaders/FlowLines.frag b/share/shaders/FlowLines.frag
index 6416bcb549..a3b7648f2e 100644
--- a/share/shaders/FlowLines.frag
+++ b/share/shaders/FlowLines.frag
@@ -57,4 +57,6 @@ void main()
         discard;
 
     fragment = color;
+    
+    if (color.a < 0.05) discard;
 }
diff --git a/share/shaders/FlowTubes.frag b/share/shaders/FlowTubes.frag
index 7afbbdab4c..0c12ffd6f3 100644
--- a/share/shaders/FlowTubes.frag
+++ b/share/shaders/FlowTubes.frag
@@ -61,5 +61,7 @@ void main() {
         discard;
 
     fragment = color;
+    
+    if (color.a < 0.05) discard;
 }