From e66f3fb04088554054e430232d0de229d93333ca Mon Sep 17 00:00:00 2001 From: Sebastian Herholz Date: Mon, 12 Oct 2026 12:28:09 +0200 Subject: [PATCH 1/5] [cleanup] Some code clean up --- openpgl/spatial/Region.h | 5 ----- tools/openpgl_debug/openpgl_debug.cpp | 1 - 2 files changed, 6 deletions(-) diff --git a/openpgl/spatial/Region.h b/openpgl/spatial/Region.h index c90e9c8b..819c1cfd 100644 --- a/openpgl/spatial/Region.h +++ b/openpgl/spatial/Region.h @@ -32,11 +32,6 @@ struct Region : public IRegion return regionBounds; } - inline const BBox &getSampleBounds() const - { - return sampleStatistics.sampleBounds; - } - #ifdef OPENPGL_RADIANCE_CACHES Vector3 getOutgoingRadiance(const Vector3 dir) const override { diff --git a/tools/openpgl_debug/openpgl_debug.cpp b/tools/openpgl_debug/openpgl_debug.cpp index 06d98abd..898009e8 100644 --- a/tools/openpgl_debug/openpgl_debug.cpp +++ b/tools/openpgl_debug/openpgl_debug.cpp @@ -36,7 +36,6 @@ struct DebugParams std::string field_file_name{""}; std::string field_file_name_comp{""}; std::string field_file_name_out{""}; - // std::string samples_file_name{""}; std::vector samples_file_names; std::string samples_file_name_comp{""}; std::string samples_out_file_name{""}; From 9f5278d2661c49a00d218b20181364334619c8f9 Mon Sep 17 00:00:00 2001 From: Sebastian Herholz Date: Mon, 12 Oct 2026 12:31:10 +0200 Subject: [PATCH 2/5] [stuff] Some minor things --- openpgl/directional/dqt/DQTFactory.h | 15 ++++++++++----- openpgl/field/Field.h | 1 + openpgl/include/openpgl/openpgl.h | 1 + 3 files changed, 12 insertions(+), 5 deletions(-) diff --git a/openpgl/directional/dqt/DQTFactory.h b/openpgl/directional/dqt/DQTFactory.h index b8837e9d..4acac5cd 100644 --- a/openpgl/directional/dqt/DQTFactory.h +++ b/openpgl/directional/dqt/DQTFactory.h @@ -164,6 +164,11 @@ class DirectionalQuadtreeFactory { return true; } + + std::string toString() const + { + return ""; + } }; struct FittingStatistics @@ -184,7 +189,7 @@ class DirectionalQuadtreeFactory const size_t numZeroValueSamples, const bool varianceBased) const {} - void fit(Distribution &dist, Statistics &stats, const SampleData *samples, const size_t numSamples, const Configuration &cfg, FittingStatistics &fitStats) + void fit(Distribution &dist, Statistics &stats, const SampleData *samples, const size_t numSamples, const Configuration &cfg, FittingStatistics &fitStats) const { for (uint32_t i = 0; i < 5; i++) { @@ -194,7 +199,7 @@ class DirectionalQuadtreeFactory } } - void update(Distribution &dist, Statistics &stats, const SampleData *samples, const size_t numSamples, const Configuration &cfg, FittingStatistics &fitStats) + void update(Distribution &dist, Statistics &stats, const SampleData *samples, const size_t numSamples, const Configuration &cfg, FittingStatistics &fitStats) const { Context ctx = {&dist, &stats, &fitStats, samples, numSamples, &cfg}; @@ -237,7 +242,7 @@ class DirectionalQuadtreeFactory // Internal Update Routines template - void updateInternal(Context &ctx) + void updateInternal(Context &ctx) const { ctx.fitStats->numSamples = ctx.stats->numSamples += ctx.numSamples; @@ -327,7 +332,7 @@ class DirectionalQuadtreeFactory } } - void buildRecursive(Context &ctx, std::vector &old_nodes, Rect rect, uint32_t i, uint32_t j, uint32_t level = 0) + void buildRecursive(Context &ctx, std::vector &old_nodes, Rect rect, uint32_t i, uint32_t j, uint32_t level = 0) const { StatsNode &old_node = old_nodes[i]; ctx.stats->nodes[j] = old_node; @@ -368,7 +373,7 @@ class DirectionalQuadtreeFactory } } - void buildSplit(Context &ctx, Rect rect, uint32_t p, uint32_t i, uint32_t level) + void buildSplit(Context &ctx, Rect rect, uint32_t p, uint32_t i, uint32_t level) const { const auto &parent = ctx.stats->nodes[p]; auto &node = ctx.stats->nodes[i]; diff --git a/openpgl/field/Field.h b/openpgl/field/Field.h index ac2020f6..a36e8aa9 100644 --- a/openpgl/field/Field.h +++ b/openpgl/field/Field.h @@ -491,6 +491,7 @@ struct Field regionStorage.first.numSplits = 0; regionStorage.first.initialized = true; } + regionStorage.first.initialized = true; } else { diff --git a/openpgl/include/openpgl/openpgl.h b/openpgl/include/openpgl/openpgl.h index 11ce5ab3..93ffd61e 100644 --- a/openpgl/include/openpgl/openpgl.h +++ b/openpgl/include/openpgl/openpgl.h @@ -18,6 +18,7 @@ #include "common.h" #include "defines.h" +#include "version.h" #if defined(OPENPGL_DIRECTION_COMPRESSION) || defined(OPENPGL_RADIANCE_COMPRESSION) #include "compression.h" #endif From 6e6cfb36b58355306042eefde1cf240f84e48617 Mon Sep 17 00:00:00 2001 From: Sebastian Herholz Date: Mon, 12 Oct 2026 12:30:37 +0200 Subject: [PATCH 3/5] [common] Adding functions to (de)serialize integer vectors --- openpgl/openpgl_common.h | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/openpgl/openpgl_common.h b/openpgl/openpgl_common.h index 1ac6423e..d1a77539 100644 --- a/openpgl/openpgl_common.h +++ b/openpgl/openpgl_common.h @@ -63,6 +63,30 @@ void deserializeFloatVectors(std::istream &stream, embree::vfloat *v } } +template +void serializeIntVectors(std::ostream &stream, const embree::vint *vectors) +{ + for (int i = 0; i < NumVectors; i++) + { + for (int j = 0; j < VectorSize; j++) + { + stream.write(reinterpret_cast(&vectors[i][j]), sizeof(int32_t)); + } + } +} + +template +void deserializeIntVectors(std::istream &stream, embree::vint *vectors) +{ + for (int i = 0; i < NumVectors; i++) + { + for (int j = 0; j < VectorSize; j++) + { + stream.read(reinterpret_cast(&vectors[i][j]), sizeof(int32_t)); + } + } +} + template void serializeVec2Vectors(std::ostream &stream, const embree::Vec2 > *vectors) { From a9480f90450ee3434fed928a1adf232385817ce0 Mon Sep 17 00:00:00 2001 From: Sebastian Herholz Date: Mon, 12 Oct 2026 12:31:45 +0200 Subject: [PATCH 4/5] [sam] Refactoring and documenting the Sam VMM fitting code --- .../vmm/AdaptiveSplitandMergeFactory.h | 85 +-- .../vmm/ParallaxAwareVonMisesFisherMixture.h | 88 ++- ...llaxAwareVonMisesFisherWeightedEMFactory.h | 560 +++++++++++------- .../vmm/VMMChiSquareComponentMerger.h | 68 ++- .../vmm/VMMChiSquareComponentSplitter.h | 15 +- 5 files changed, 528 insertions(+), 288 deletions(-) diff --git a/openpgl/directional/vmm/AdaptiveSplitandMergeFactory.h b/openpgl/directional/vmm/AdaptiveSplitandMergeFactory.h index c676ff62..c32f91b7 100644 --- a/openpgl/directional/vmm/AdaptiveSplitandMergeFactory.h +++ b/openpgl/directional/vmm/AdaptiveSplitandMergeFactory.h @@ -29,13 +29,15 @@ struct AdaptiveSplitAndMergeFactory // typedef WeightedEMVonMisesFisherFactory WeightedEMFactory; typedef ParallaxAwareVonMisesFisherWeightedEMFactory WeightedEMFactory; typedef VonMisesFisherChiSquareComponentSplitter Splitter; - typedef VonMisesFisherChiSquareComponentMerger Merger; + typedef VonMisesFisherChiSquareComponentMerger Merger; struct Configuration { typename WeightedEMFactory::Configuration weightedEMCfg; + // The min. Chi^2 threshold for splitting a mixture component float splittingThreshold{0.5f}; + // The max. Chi^2 threshold for merging a two mixture components float mergingThreshold{0.025f}; bool useSplitAndMerge{true}; @@ -43,8 +45,11 @@ struct AdaptiveSplitAndMergeFactory bool partialReFit{false}; int maxSplitItr{1}; + // The min. number of samples processed after the last merge step required to trigger a splitting step int minSamplesForSplitting{0}; + // The min. number of samples processed after the last merge step required to trigger partial refitting after splitting int minSamplesForPartialRefitting{0}; + // The min. number of samples processed after the last merge step required to trigger a merge step int minSamplesForMerging{0}; void serialize(std::ostream &stream) const; @@ -71,8 +76,6 @@ struct AdaptiveSplitAndMergeFactory typename WeightedEMFactory::SufficientStatistics sufficientStatistics; typename Splitter::ComponentSplitStatistics splittingStatistics; - // size_t numComponents {0}; - size_t numSamplesAfterLastSplit{0}; size_t numSamplesAfterLastMerge{0}; @@ -308,7 +311,7 @@ void AdaptiveSplitAndMergeFactory::fit(VMM &vmm, Statistics &s { const size_t numComponents = cfg.weightedEMCfg.initK; stats.clear(numComponents); - // intial fit + // Initial fitting of the mixture using standard weighted EM WeightedEMFactory factory = WeightedEMFactory(); typename WeightedEMFactory::FittingStatistics wemFitStats; factory.fitMixture(vmm, stats.sufficientStatistics, samples, numSamples, cfg.weightedEMCfg, wemFitStats); @@ -316,28 +319,29 @@ void AdaptiveSplitAndMergeFactory::fit(VMM &vmm, Statistics &s OPENPGL_ASSERT(vmm.isValid()); OPENPGL_ASSERT(vmm.getNumComponents() == stats.sufficientStatistics.getNumComponents()); OPENPGL_ASSERT(stats.isValid()); - /* */ + // We use split and merge to optimze the fitting result of the standard weighted EM algorithm which can get stuck in local + // maximas (e.g., one component trying to represent a multi-modal distribution or a distribtion containing a firely signal). if (cfg.useSplitAndMerge) { - // calculate the estimate of the integral of the function (e.g. radiance or importance) fitted by the VMM + // Calculate the estimate of the integral of the function (e.g. radiance or importance) represented by the VMM float mcEstimate = stats.sufficientStatistics.getSumWeights() / stats.sufficientStatistics.getNumSamples(); - // split the fitted components of the inital fit to match - // the observed samples -#ifdef OPENPGL_SHOW_PRINT_OUTS - std::cout << stats.sufficientStatistics.toString() << std::endl; -#endif + // Perfroming recursive splitting so that the mixture components best represent the observed data Splitter splitter = Splitter(); splitter.PerformRecursiveSplitting(vmm, stats.sufficientStatistics, cfg.splittingThreshold, mcEstimate, samples, numSamples, cfg.weightedEMCfg); + // After the recursive splitting process is finished we recaclulate the split statistics for current mixture from scratch splitter.CalculateSplitStatistics(vmm, stats.splittingStatistics, mcEstimate, samples, numSamples); OPENPGL_ASSERT(vmm.getNumComponents() == stats.getNumComponents()); OPENPGL_ASSERT(vmm.isValid()); + ////////////////////////////////////////////////////// + // Merging + ////////////////////////////////////////////////////// Merger merger = Merger(); - merger.PerformMerging(vmm, cfg.mergingThreshold, stats.sufficientStatistics, stats.splittingStatistics); + merger.PerformMerging(vmm, cfg.mergingThreshold, cfg.splittingThreshold, false, stats.sufficientStatistics, stats.splittingStatistics); OPENPGL_ASSERT(vmm.isValid()); } @@ -357,24 +361,25 @@ void AdaptiveSplitAndMergeFactory::update(VMM &vmm, Statistics OPENPGL_ASSERT(vmm.getNumComponents() == stats.getNumComponents()); OPENPGL_ASSERT(stats.isValid()); - // first update the mixture + // Update the mixture using standard weighted EM WeightedEMFactory factory = WeightedEMFactory(); typename WeightedEMFactory::FittingStatistics wemFitStats; - // stats.sufficientStatistics.clear(vmm._numComponents); const size_t prevNumberOfComponents = vmm._numComponents; factory.updateMixture(vmm, stats.sufficientStatistics, samples, numSamples, cfg.weightedEMCfg, wemFitStats); OPENPGL_ASSERT(vmm.isValid()); - // check if the update step added a new component. - // This happnes if samples are not covered by any existing component - // and we need to extend the splittingStats. + + // Check if the update step added a new component. + // This happens if samples are not covered by any existing component and we need to extend the splittingStats. if (prevNumberOfComponents < vmm._numComponents) { stats.splittingStatistics.setNumComponents(vmm._numComponents); } OPENPGL_ASSERT(stats.sufficientStatistics.isValid()); + // We use split and merge to optimize the fitting/update of the mixture to better reflect the observed data. if (cfg.useSplitAndMerge) { + // Calculate the estimate of the integral of the function (e.g. radiance or importance) represented by the VMM float mcEstimate = stats.sufficientStatistics.getSumWeights() / stats.sufficientStatistics.getNumSamples(); fitStats.numSamples = numSamples; @@ -383,28 +388,42 @@ void AdaptiveSplitAndMergeFactory::update(VMM &vmm, Statistics stats.numSamplesAfterLastSplit += numSamples; stats.numSamplesAfterLastMerge += numSamples; + ////////////////////////////////////////////////////// + // Splitting + ////////////////////////////////////////////////////// Splitter splitter = Splitter(); - // OPENPGL_ASSERT(stats.splittingStatistics.isValid()); + // Updating split statistics splitter.UpdateSplitStatistics(vmm, stats.splittingStatistics, mcEstimate, samples, numSamples); OPENPGL_ASSERT(stats.splittingStatistics.isValid()); + // We only perform splitting if we have observed enough samples after the last splitting round. if (stats.numSamplesAfterLastSplit >= cfg.minSamplesForSplitting) { + // The binary mask tagging the split components which need to be refitted typename WeightedEMFactory::PartialFittingMask mask; mask.resetToFalse(); + // The binary mask identifying if the previous components stats of the split components should be used as prior or not. + // In this version it is alsways set to false. + typename WeightedEMFactory::PartialFittingMask previousAsPriorMask; + previousAsPriorMask.resetToFalse(); + // Getting the list of split candidates sorted by their chi^2 values std::vector splitComps = stats.splittingStatistics.getSplitCandidates(); int totalSplitCount = 0; - // const size_t numComp = vmm._numComponents; + // For each split cadidate we check if its chi^2 value is above our split threshold and if we still have free components in our mixture. for (size_t k = 0; k < splitComps.size(); k++) { if (splitComps[k].chiSquareEst > cfg.splittingThreshold && vmm._numComponents < VMM::MaxComponents) { - splitter.SplitComponent(vmm, stats.splittingStatistics, stats.sufficientStatistics, splitComps[k].componentIndex); - mask.setToTrue(splitComps[k].componentIndex); - mask.setToTrue(vmm._numComponents - 1); - // std::cout << "split[" << totalSplitCount << "]: " << "\tidx0: " << splitComps[k].componentIndex << "\tidx1: " << vmm._numComponents-1 << std::endl; - totalSplitCount++; + // We split the component. + bool success = splitter.SplitComponent(vmm, stats.splittingStatistics, stats.sufficientStatistics, splitComps[k].componentIndex); + // If the split was successfull mark the new components for partial refitting. + if (success) + { + mask.setToTrue(splitComps[k].componentIndex); + mask.setToTrue(vmm._numComponents - 1); + totalSplitCount++; + } } } @@ -412,15 +431,12 @@ void AdaptiveSplitAndMergeFactory::update(VMM &vmm, Statistics OPENPGL_ASSERT(vmm.getNumComponents() == stats.getNumComponents()); OPENPGL_ASSERT(stats.isValid()); + // We perform a partial refit, if we performed any split and the number of samples is above our minimal sample threshold if (totalSplitCount > 0 && cfg.partialReFit && numSamples >= cfg.minSamplesForPartialRefitting) { - typename WeightedEMFactory::SufficientStatistics tempSuffStatistics = stats.sufficientStatistics; - tempSuffStatistics.clear(vmm._numComponents); - factory.partialUpdateMixture(vmm, mask, tempSuffStatistics, samples, numSamples, cfg.weightedEMCfg, wemFitStats); - stats.sufficientStatistics.setNumComponents(vmm._numComponents); - stats.sufficientStatistics.maskedReplace(mask, tempSuffStatistics); - // update number of components for the splitStats to - // account for additionaly added componetes based on not covered samples. + // For this partial refitting we are NOT using the previous component stats as prior + factory.partialUpdateMixture(vmm, mask, false, previousAsPriorMask, stats.sufficientStatistics, samples, numSamples, cfg.weightedEMCfg, wemFitStats); + // Updating the number of component in the splitting stats to account for the new split components stats.splittingStatistics.setNumComponents(vmm._numComponents); fitStats.numPartialUpdateWEMIterations = wemFitStats.numIterations; OPENPGL_ASSERT(vmm.isValid()); @@ -435,15 +451,18 @@ void AdaptiveSplitAndMergeFactory::update(VMM &vmm, Statistics std::cout << "update: totalSplitCount = " << totalSplitCount << "\t splitThreshold: " << cfg.splittingThreshold << std::endl; #endif } - + // OLD OPENPGL_ASSERT(vmm.isValid()); OPENPGL_ASSERT(vmm.getNumComponents() == stats.getNumComponents()); OPENPGL_ASSERT(stats.isValid()); + ////////////////////////////////////////////////////// + // Merging + ////////////////////////////////////////////////////// if (stats.numSamplesAfterLastMerge >= cfg.minSamplesForMerging) { Merger merger = Merger(); - size_t numMerges = merger.PerformMerging(vmm, cfg.mergingThreshold, stats.sufficientStatistics, stats.splittingStatistics); + size_t numMerges = merger.PerformMerging(vmm, cfg.mergingThreshold, cfg.splittingThreshold, false, stats.sufficientStatistics, stats.splittingStatistics); fitStats.numMerges = numMerges; stats.numSamplesAfterLastMerge = 0.0f; diff --git a/openpgl/directional/vmm/ParallaxAwareVonMisesFisherMixture.h b/openpgl/directional/vmm/ParallaxAwareVonMisesFisherMixture.h index 55951fc2..b23bfa33 100644 --- a/openpgl/directional/vmm/ParallaxAwareVonMisesFisherMixture.h +++ b/openpgl/directional/vmm/ParallaxAwareVonMisesFisherMixture.h @@ -99,6 +99,8 @@ struct ParallaxAwareVonMisesFisherMixture float pdf(Vector3 direction) const; + float pdfComponent(Vector3 direction, const size_t idx) const; + Vector3 sample(const Vector2 sample) const; #ifdef USE_SIMD_CDF_SAMPLING @@ -111,6 +113,8 @@ struct ParallaxAwareVonMisesFisherMixture void splitComponent(const size_t &idx0, const size_t &idx1, const float &weight0, const float &weight1, const Vector3 &meanDirection0, const Vector3 &meanDirection1, const float &meanCosine0, const float &meanCosine1); + void splitFireFlyComponent(const size_t &idx0, const size_t &idx1, const float &fireFlyFrac, const Vector3 &fireFlyMeanDirection, const float &fireFlyMeanCosine); + void performRelativeParallaxShift(const Vector3 &shiftDirection); #ifdef OPENPGL_RADIANCE_CACHES @@ -135,6 +139,8 @@ struct ParallaxAwareVonMisesFisherMixture void clearComponent(const size_t &idx); + void clearComponents(); + // Getter methods for the PAVMM attributes size_t getNumComponents() const; @@ -245,8 +251,8 @@ std::string ParallaxAwareVonMisesFisherMixture_numComponents << std::endl; float sumWeights = 0.0f; - // for ( int k = 0; k < this->_numComponents; k++) - for (int k = 0; k < maxComponents; k++) + for (int k = 0; k < this->_numComponents; k++) + // for (int k = 0; k < maxComponents; k++) { const div_t tmp = div(k, static_cast(VecSize)); ss << "vmm[" << k << "]: " << "weight: " << this->_weights[tmp.quot][tmp.rem]; @@ -343,6 +349,30 @@ void ParallaxAwareVonMisesFisherMixture +void ParallaxAwareVonMisesFisherMixture::splitFireFlyComponent(const size_t &idx0, const size_t &idx1, const float &fireFlyFrac, + const Vector3 &fireFlyMeanDirection, const float &fireFlyMeanCosine) +{ + const div_t tmpIdx0 = div(idx0, static_cast(VectorSize)); + const div_t tmpIdx1 = div(idx1, static_cast(VectorSize)); + + const float weight = _weights[tmpIdx0.quot][tmpIdx0.rem]; + _weights[tmpIdx0.quot][tmpIdx0.rem] *= (1.f - fireFlyFrac); + + _weights[tmpIdx1.quot][tmpIdx1.rem] = weight * fireFlyFrac; + _meanDirections[tmpIdx1.quot].x[tmpIdx1.rem] = fireFlyMeanDirection.x; + _meanDirections[tmpIdx1.quot].y[tmpIdx1.rem] = fireFlyMeanDirection.y; + _meanDirections[tmpIdx1.quot].z[tmpIdx1.rem] = fireFlyMeanDirection.z; + _meanCosines[tmpIdx1.quot][tmpIdx1.rem] = fireFlyMeanCosine; + _kappas[tmpIdx1.quot][tmpIdx1.rem] = MeanCosineToKappa(fireFlyMeanCosine); + _distances[tmpIdx1.quot][tmpIdx1.rem] = _distances[tmpIdx0.quot][tmpIdx0.rem]; + if (idx1 == _numComponents) + { + _numComponents++; + } + _calculateNormalization(); +} + template void ParallaxAwareVonMisesFisherMixture::mergeComponents(const size_t &idx0, const size_t &idx1) { @@ -517,6 +547,39 @@ void ParallaxAwareVonMisesFisherMixture +void ParallaxAwareVonMisesFisherMixture::clearComponents() +{ + const embree::vfloat ones(1.0f); + const embree::vfloat zeros(0.0f); + const embree::vfloat zeroKappaNorm(ONE_OVER_FOUR_PI); + + for (uint32_t k = 0; k < NumVectors; k++) + { + _weights[k] = zeros; + _kappas[k] = zeros; + _eMinus2Kappa[k] = ones; + _meanCosines[k] = zeros; + _normalizations[k] = zeroKappaNorm; + + _meanDirections[k].x = zeros; + _meanDirections[k].y = zeros; + _meanDirections[k].z = ones; + + _distances[k] = zeros; + +#ifdef OPENPGL_RADIANCE_CACHES + _fluenceRGBWeightsWithMIS[k].x = zeros; + _fluenceRGBWeightsWithMIS[k].y = zeros; + _fluenceRGBWeightsWithMIS[k].z = zeros; + + _fluenceRGBWeights[k].x = zeros; + _fluenceRGBWeights[k].y = zeros; + _fluenceRGBWeights[k].z = zeros; +#endif + } +} + template void ParallaxAwareVonMisesFisherMixture::serialize(std::ostream &stream) const { @@ -879,6 +942,27 @@ float ParallaxAwareVonMisesFisherMixture +float ParallaxAwareVonMisesFisherMixture::pdfComponent(Vector3 direction, const size_t idx) const +{ + const div_t tmpIdx = div(idx, VecSize); + + embree::vfloat pdf = {0.0f}; + embree::Vec3> vec3Direction(direction[0], direction[1], direction[2]); + + const embree::vfloat ones(1.0f); + const embree::vfloat zeros(0.0f); + + size_t k = tmpIdx.quot; + + const embree::vfloat cosTheta = embree::dot(vec3Direction, _meanDirections[k]); + const embree::vfloat cosThetaMinusOne = embree::min(cosTheta - ones, zeros); + const embree::vfloat eval = _normalizations[k] * embree::fastapprox::exp>(_kappas[k] * cosThetaMinusOne); + pdf += _weights[k] * eval; + + return pdf[tmpIdx.rem]; +} + template bool ParallaxAwareVonMisesFisherMixture::softAssignment( Vector3 direction, typename ParallaxAwareVonMisesFisherMixture::SoftAssignment &softAssign) const diff --git a/openpgl/directional/vmm/ParallaxAwareVonMisesFisherWeightedEMFactory.h b/openpgl/directional/vmm/ParallaxAwareVonMisesFisherWeightedEMFactory.h index 033c67b2..1e85d95d 100644 --- a/openpgl/directional/vmm/ParallaxAwareVonMisesFisherWeightedEMFactory.h +++ b/openpgl/directional/vmm/ParallaxAwareVonMisesFisherWeightedEMFactory.h @@ -83,6 +83,7 @@ struct ParallaxAwareVonMisesFisherWeightedEMFactory PartialFittingMask() = default; void resetToFalse(); + void resetToTrue(); void resetToTrue(const size_t &numComponents); void setToTrue(const size_t &idx); void setToFalse(const size_t &idx); @@ -90,7 +91,7 @@ struct ParallaxAwareVonMisesFisherWeightedEMFactory std::string toString() const; }; - struct SufficientStatistics //: public WEMVMMFactory::SufficientStatistics + struct SufficientStatistics { // FittingStatistics public: @@ -103,11 +104,22 @@ struct ParallaxAwareVonMisesFisherWeightedEMFactory embree::vfloat surfaceSampleNumberWeights[VMM::NumVectors]; #endif float sumWeights{0.f}; + + // number of effective samples for the fit this number can be smaller than overallNumSamples + // due to decaying the statistics float numSamples{0.f}; + + // the total number of samples used for fitting without the decay float overallNumSamples{0.f}; + + // Number of mixture components size_t numComponents{VMM::MaxComponents}; + // If the statistics are already normalized or not bool normalized{false}; + float norm{1.f}; + float inv_norm{1.f}; + embree::vfloat sumOfDistanceWeightes[VMM::NumVectors]; SufficientStatistics() = default; @@ -136,6 +148,8 @@ struct ParallaxAwareVonMisesFisherWeightedEMFactory void splitComponentsStats(const size_t &idx0, const size_t &idx1, const Vector3 &meanDirection0, const Vector3 &meanDirection1, const float &meanCosine0, const float &meanCosine1); + void splitFireFlyComponentsStats(const size_t &idx0, const size_t &idx1, const float fireFlyFrac, const Vector3 &fireFlyMeanDirection, const float &fireFlyMeanCosine); + void swapComponentStats(const size_t &idx0, const size_t &idx1); void maskedReplace(const PartialFittingMask &mask, const SufficientStatistics &stats); @@ -154,6 +168,11 @@ struct ParallaxAwareVonMisesFisherWeightedEMFactory return sumWeights; } + inline float getMeanSamplesWeights() const + { + return sumWeights / numSamples; + } + inline void setNumComponents(const size_t &numComponents) { this->numComponents = numComponents; @@ -191,8 +210,11 @@ struct ParallaxAwareVonMisesFisherWeightedEMFactory void updateMixture(VMM &vmm, SufficientStatistics &previousStats, const SampleData *samples, const size_t numSamples, const Configuration &cfg, FittingStatistics &fitStats) const; - void partialUpdateMixture(VMM &vmm, PartialFittingMask &mask, SufficientStatistics &previousStats, const SampleData *samples, const size_t numSamples, const Configuration &cfg, - FittingStatistics &fitStats) const; + void partialUpdateMixture(VMM &vmm, PartialFittingMask &mask, const bool usePreviousStatsAsPrior, PartialFittingMask &previousAsPriorMask, SufficientStatistics &previousStats, + const SampleData *samples, const size_t numSamples, const Configuration &cfg, FittingStatistics &fitStats) const; + + void partialMergeSufficientStatisticsWithPriors(PartialFittingMask &mask, SufficientStatistics ¤tStats, const bool usePreviousStatsAsPrior, + const PartialFittingMask &previousAsPriorMask, const SufficientStatistics &previousStats) const; #ifdef OPENPGL_RADIANCE_CACHES void updateFluenceEstimate(VMM &vmm, const SampleData *samples, const size_t numSamples, const size_t numZeroValueSamples, const SampleStatistics &sampleStatistics) const; @@ -218,20 +240,11 @@ struct ParallaxAwareVonMisesFisherWeightedEMFactory float weightedExpectationStep(VMM &vmm, SufficientStatistics &stats, UnassignedSamplesStatistics &unassignedStats, const SampleData *samples, const size_t numSamples) const; - void weightedMaximumAPosteriorStep(VMM &vmm, const SufficientStatistics &previousStats, const SufficientStatistics ¤tStats, const Configuration &cfg) const; - - void estimateMAPWeights(VMM &vmm, const SufficientStatistics ¤tStats, const SufficientStatistics &previousStats, const float &_weightPrior) const; - - void estimateMAPMeanDirectionAndConcentration(VMM &vmm, const SufficientStatistics ¤tStats, const SufficientStatistics &previousStats, const Configuration &cfg) const; - - void partialWeightedMaximumAPosteriorStep(VMM &vmm, const PartialFittingMask &mask, SufficientStatistics &previousStats, SufficientStatistics ¤tStats, - const Configuration &cfg) const; + void weightedMaximumAPosteriorStep(VMM &vmm, const SufficientStatistics ¤tStats, const Configuration &cfg) const; - void estimatePartialMAPWeights(VMM &vmm, const PartialFittingMask &mask, SufficientStatistics ¤tStats, SufficientStatistics &previousStats, - const float &_weightPrior) const; + void estimateMAPWeights(VMM &vmm, const SufficientStatistics ¤tStats, const float &_weightPrior) const; - void estimatePartialMAPMeanDirectionAndConcentration(VMM &vmm, const PartialFittingMask &mask, SufficientStatistics ¤tStats, SufficientStatistics &previousStats, - const Configuration &cfg) const; + void estimateMAPMeanDirectionAndConcentration(VMM &vmm, const SufficientStatistics ¤tStats, const Configuration &cfg) const; void handleUnassignedSampleStats(UnassignedSamplesStatistics &unassignedStats, VMM &vmm, SufficientStatistics ¤tStats, SufficientStatistics &previousStats) const; @@ -403,6 +416,8 @@ void ParallaxAwareVonMisesFisherWeightedEMFactory::SufficientS stream.write(reinterpret_cast(&overallNumSamples), sizeof(float)); stream.write(reinterpret_cast(&numComponents), sizeof(size_t)); stream.write(reinterpret_cast(&normalized), sizeof(bool)); + stream.write(reinterpret_cast(&norm), sizeof(float)); + stream.write(reinterpret_cast(&inv_norm), sizeof(float)); #ifdef OPENPGL_VSP_GUIDING serializeFloatVectors(stream, volumeContributionWeights); serializeFloatVectors(stream, surfaceContributionWeights); @@ -422,6 +437,8 @@ void ParallaxAwareVonMisesFisherWeightedEMFactory::SufficientS stream.read(reinterpret_cast(&overallNumSamples), sizeof(float)); stream.read(reinterpret_cast(&numComponents), sizeof(size_t)); stream.read(reinterpret_cast(&normalized), sizeof(bool)); + stream.read(reinterpret_cast(&norm), sizeof(float)); + stream.read(reinterpret_cast(&inv_norm), sizeof(float)); #ifdef OPENPGL_VSP_GUIDING deserializeFloatVectors(stream, volumeContributionWeights); deserializeFloatVectors(stream, surfaceContributionWeights); @@ -454,8 +471,11 @@ void ParallaxAwareVonMisesFisherWeightedEMFactory::SufficientS } sumWeights = 0.0f; - numSamples = 0; + numSamples = 0.0f; + overallNumSamples = 0.0f; normalized = false; + norm = 1.0f; + inv_norm = 1.0f; } template @@ -510,6 +530,7 @@ template std::string ParallaxAwareVonMisesFisherWeightedEMFactory::SufficientStatistics::toString() const { std::stringstream ss; + ss << std::setprecision(10); ss << "SufficientStatistics:" << std::endl; ss << "\tsumWeights = " << sumWeights << std::endl; ss << "\tnumSamples = " << numSamples << std::endl; @@ -517,13 +538,21 @@ std::string ParallaxAwareVonMisesFisherWeightedEMFactory::Suff ss << "\tnumComponents = " << numComponents << std::endl; ss << "\tisNormalized = " << normalized << std::endl; // for (size_t k = 0; k < numComponents ; k++) + + ss << "\tnorm = " << norm << std::endl; + ss << "\tinv_norm = " << inv_norm << std::endl; + + float sumWeightedStats{0.0f}; for (size_t k = 0; k < VMM::MaxComponents; k++) { int i = k / VMM::VectorSize; int j = k % VMM::VectorSize; - ss << "\tstat[" << k << "]:" << "\tsumWeightedStats = " << sumOfWeightedStats[i][j] << "\tsumWeightedDirections = [" << sumOfWeightedDirections[i].x[j] << ",\t" - << sumOfWeightedDirections[i].y[j] << ",\t" << sumOfWeightedDirections[i].z[j] << "]" << "\tsumWeightedDistanceWeights = " << sumOfDistanceWeightes[i][j] << std::endl; + ss << "\tstat[" << k << "]:" << "\tsumWeightedStats = " << sumOfWeightedStats[i][j] << "\tsumWeightedStats = " << sumOfWeightedStats[i][j] * inv_norm + << "\tsumWeightedDirections = [" << sumOfWeightedDirections[i].x[j] << ",\t" << sumOfWeightedDirections[i].y[j] << ",\t" << sumOfWeightedDirections[i].z[j] << "]" + << "\tsumWeightedDistanceWeights = " << sumOfDistanceWeightes[i][j] << std::endl; + sumWeightedStats += sumOfWeightedStats[i][j]; } + ss << "\tsumWeightedStats = " << sumWeightedStats << std::endl; return ss.str(); } @@ -548,6 +577,9 @@ void ParallaxAwareVonMisesFisherWeightedEMFactory::SufficientS { sumWeights = embree::reduce_add(newSumWeights); } + + norm = numSamples / sumWeights; + inv_norm = sumWeights / numSamples; } template @@ -643,6 +675,35 @@ void ParallaxAwareVonMisesFisherWeightedEMFactory::SufficientS OPENPGL_ASSERT(!std::isnan(sumOfWeightedDirections[tmpI.quot].z[tmpI.rem]) && std::isfinite(sumOfWeightedDirections[tmpI.quot].z[tmpI.rem])); } +template +void ParallaxAwareVonMisesFisherWeightedEMFactory::SufficientStatistics::splitFireFlyComponentsStats(const size_t &idx0, const size_t &idx1, + const float fireFlyFrac, const Vector3 &fireFlyMeanDirection, + const float &fireFlyMeanCosine) +{ + const div_t tmpI = div(idx0, static_cast(VMM::VectorSize)); + const div_t tmpJ = div(idx1, static_cast(VMM::VectorSize)); + + sumOfWeightedStats[tmpJ.quot][tmpJ.rem] = sumOfWeightedStats[tmpI.quot][tmpI.rem]; + sumOfWeightedStats[tmpI.quot][tmpI.rem] *= 1.f - fireFlyFrac; + sumOfWeightedStats[tmpJ.quot][tmpJ.rem] *= fireFlyFrac; + + sumOfWeightedDirections[tmpI.quot].x[tmpI.rem] *= 1.f - fireFlyFrac; + sumOfWeightedDirections[tmpI.quot].y[tmpI.rem] *= 1.f - fireFlyFrac; + sumOfWeightedDirections[tmpI.quot].z[tmpI.rem] *= 1.f - fireFlyFrac; + + Vector3 newMeanDirection = fireFlyMeanDirection; + newMeanDirection *= fireFlyMeanCosine * sumOfWeightedStats[tmpJ.quot][tmpJ.rem]; + sumOfWeightedDirections[tmpJ.quot].x[tmpJ.rem] = newMeanDirection.x; + sumOfWeightedDirections[tmpJ.quot].y[tmpJ.rem] = newMeanDirection.y; + sumOfWeightedDirections[tmpJ.quot].z[tmpJ.rem] = newMeanDirection.z; + + numComponents += 1; + + OPENPGL_ASSERT(!std::isnan(sumOfWeightedDirections[tmpI.quot].x[tmpI.rem]) && std::isfinite(sumOfWeightedDirections[tmpI.quot].x[tmpI.rem])); + OPENPGL_ASSERT(!std::isnan(sumOfWeightedDirections[tmpI.quot].y[tmpI.rem]) && std::isfinite(sumOfWeightedDirections[tmpI.quot].y[tmpI.rem])); + OPENPGL_ASSERT(!std::isnan(sumOfWeightedDirections[tmpI.quot].z[tmpI.rem]) && std::isfinite(sumOfWeightedDirections[tmpI.quot].z[tmpI.rem])); +} + template void ParallaxAwareVonMisesFisherWeightedEMFactory::SufficientStatistics::normalize(const float &_numSamples) { @@ -655,35 +716,43 @@ void ParallaxAwareVonMisesFisherWeightedEMFactory::SufficientS sumWeightedStatsVec += sumOfWeightedStats[k]; } sumWeights = reduce_add(sumWeightedStatsVec); - embree::vfloat norm(sumWeights > FLT_EPSILON ? _numSamples / sumWeights : 1.0f); + norm = sumWeights > FLT_EPSILON ? _numSamples / sumWeights : 1.0f; + inv_norm = _numSamples > FLT_EPSILON ? sumWeights / _numSamples : 1.0f; + embree::vfloat normVec(_numSamples / sumWeights); for (int k = 0; k < cnt; k++) { - sumOfWeightedDirections[k] *= norm; - sumOfWeightedStats[k] *= norm; + sumOfWeightedDirections[k] *= normVec; + sumOfWeightedStats[k] *= normVec; } normalized = true; + // TODO: check this seems to be wrong + overallNumSamples = numSamples; } - +/* */ template typename ParallaxAwareVonMisesFisherWeightedEMFactory::SufficientStatistics & ParallaxAwareVonMisesFisherWeightedEMFactory::SufficientStatistics::operator+=( const typename ParallaxAwareVonMisesFisherWeightedEMFactory::SufficientStatistics &stats) { + OPENPGL_ASSERT(this->numComponents == stats.numComponents); // TODO: check for normalization + if ((this->overallNumSamples > 0.f && !this->normalized) || (stats.overallNumSamples > 0.f && !stats.normalized)) + std::cout << "ERROR: normalization" << std::endl; const int cnt = (numComponents + VMM::VectorSize - 1) / VMM::VectorSize; this->sumWeights += stats.sumWeights; this->numSamples += stats.numSamples; - this->overallNumSamples += stats.numSamples; + this->overallNumSamples += stats.overallNumSamples; for (int k = 0; k < cnt; k++) { this->sumOfWeightedDirections[k] += stats.sumOfWeightedDirections[k]; this->sumOfWeightedStats[k] += stats.sumOfWeightedStats[k]; this->sumOfDistanceWeightes[k] += stats.sumOfDistanceWeightes[k]; } - + this->norm = this->numSamples / this->sumWeights; + this->inv_norm = this->sumWeights / this->numSamples; return *this; } @@ -793,11 +862,10 @@ template typename ParallaxAwareVonMisesFisherWeightedEMFactory::VMM ParallaxAwareVonMisesFisherWeightedEMFactory::VMMfromSufficientStatistics( const SufficientStatistics &suffStats, const Configuration &cfg) const { - SufficientStatistics previousStats; - previousStats.clear(suffStats.numComponents); VMM vmm; + // vmm.clearComponents(); // TODO: seems to be a debug check vmm._numComponents = suffStats.numComponents; - weightedMaximumAPosteriorStep(vmm, previousStats, suffStats, cfg); + weightedMaximumAPosteriorStep(vmm, suffStats, cfg); return vmm; } @@ -807,13 +875,9 @@ void ParallaxAwareVonMisesFisherWeightedEMFactory::fitMixture( const Configuration &cfg, FittingStatistics &fitStats) const { const size_t numComponents = cfg.initK; - // VonMisesFisherFactory< TVMMDistribution>::InitUniformVMM( vmm, numComponents, 5.0f); this->InitUniformVMM(vmm, numComponents, cfg.initKappa); - // OPENPGL_ASSERT(vmm.isValid()); - // SufficientStatistics stats; stats.clear(numComponents); stats.normalized = true; - // stats.clearAll(); updateMixture(vmm, stats, samples, numSamples, cfg, fitStats); } @@ -821,6 +885,11 @@ template void ParallaxAwareVonMisesFisherWeightedEMFactory::handleUnassignedSampleStats(UnassignedSamplesStatistics &unassignedStats, VMM &vmm, SufficientStatistics ¤tStats, SufficientStatistics &previousStats) const { + OPENPGL_ASSERT(embree::isvalid(unassignedStats.sumOfUnassignedWeights)); + OPENPGL_ASSERT(embree::isvalid(unassignedStats.sumUnassignedWeightedDirections.x)); + OPENPGL_ASSERT(embree::isvalid(unassignedStats.sumUnassignedWeightedDirections.y)); + OPENPGL_ASSERT(embree::isvalid(unassignedStats.sumUnassignedWeightedDirections.z)); + const div_t tmpK = div(currentStats.numComponents, TVMMDistribution::VectorSize); currentStats.numComponents++; currentStats.sumOfWeightedStats[tmpK.quot][tmpK.rem] = unassignedStats.sumOfUnassignedWeights; @@ -850,21 +919,30 @@ void ParallaxAwareVonMisesFisherWeightedEMFactory::updateMixtu float previousLogLikelihood = 0.0f; float inv_previousLogLikelihood = 1.0f; UnassignedSamplesStatistics unassignedStats; + + // Running multiple EM iterations until the mixture is converged or a number of max iterations is reached while (!converged && currentEMIteration < cfg.maxEMIterrations) { + // Running the E-step to calculate the sufficient statistics and estimate the current log likelihood float logLikelihood = weightedExpectationStep(vmm, currentStats, unassignedStats, samples, numSamples); + // Special handling of samples which are not covered by any mixture component (i.e., adding an additional/special component) if (unassignedStats.sumOfUnassignedWeights > 0.0f && currentStats.numComponents < TVMMDistribution::MaxComponents) { handleUnassignedSampleStats(unassignedStats, vmm, currentStats, previousStats); } OPENPGL_ASSERT(!currentStats.isNormalized()); + // Normalizing sufficient statistics so that the weighted stats per component can be re-interpreded by number of samples currentStats.normalize(currentStats.numSamples); OPENPGL_ASSERT(currentStats.isValid()); - weightedMaximumAPosteriorStep(vmm, currentStats, previousStats, cfg); + // Adding/Merging the sufficient statistics of the previous training/update iteration to act as a prior to implement + // incremental EM fitting + currentStats += previousStats; + // Deriving the mixture parameter from the merged suffcient statistics using MAP and addtional global priors. + weightedMaximumAPosteriorStep(vmm, currentStats, cfg); currentEMIteration++; - // TODO: Add convergence check + // Checking of the fitting converged (i.e., the loglikelihood only changed minimally) if (currentEMIteration > 1) { float relLogLikelihoodDifference = std::fabs(logLikelihood - previousLogLikelihood) * inv_previousLogLikelihood; @@ -872,22 +950,22 @@ void ParallaxAwareVonMisesFisherWeightedEMFactory::updateMixtu { converged = true; } - // std::cout << "logLikelihood:" << logLikelihood << "\t previousLogLikelihood: "<< previousLogLikelihood << "\t relLogLikelihoodDifference: " << - // relLogLikelihoodDifference << std::endl; previousLogLikelihood = logLikelihood; inv_previousLogLikelihood = 1.0f / std::fabs(logLikelihood); } } - previousStats += currentStats; + + // The merged sufficient stats from the last iteration are now the new previous/prior stats + previousStats = currentStats; fitStats.numSamples = numSamples; fitStats.numIterations = currentEMIteration; fitStats.summedWeightedLogLikelihood = previousLogLikelihood; - // td::cout << "converged:" << currentEMIteration << std::endl; } template -void ParallaxAwareVonMisesFisherWeightedEMFactory::partialUpdateMixture(VMM &vmm, PartialFittingMask &mask, SufficientStatistics &previousStats, +void ParallaxAwareVonMisesFisherWeightedEMFactory::partialUpdateMixture(VMM &vmm, PartialFittingMask &mask, const bool usePreviousStatsAsPrior, + PartialFittingMask &previousAsPriorMask, SufficientStatistics &previousStats, const SampleData *samples, const size_t numSamples, const Configuration &cfg, FittingStatistics &fitStats) const { @@ -899,25 +977,33 @@ void ParallaxAwareVonMisesFisherWeightedEMFactory::partialUpda bool converged = false; float previousLogLikelihood = 0.0f; float inv_previousLogLikelihood = 1.0f; - UnassignedSamplesStatistics unassignedStats; + + // Running multiple EM iterations until the mixture is converged or a number of max iterations is reached. + // During these iterations only the masked mixture components are updated while (!converged && currentEMIteration < cfg.maxEMIterrations) { + // Running the E-step to calculate the sufficient statistics and estimate the current log likelihood float logLikelihood = weightedExpectationStep(vmm, currentStats, unassignedStats, samples, numSamples); + // Special handling of samples which are not covered by any mixture component (i.e., adding an additional/special component) if (unassignedStats.sumOfUnassignedWeights > 0.0f && currentStats.numComponents < TVMMDistribution::MaxComponents) { handleUnassignedSampleStats(unassignedStats, vmm, currentStats, previousStats); mask.setToTrue(vmm._numComponents - 1); } - OPENPGL_ASSERT(currentStats.isValid()); OPENPGL_ASSERT(!currentStats.isNormalized()); + // Normalizing sufficient statistics so that the weighted stats per component can be re-interpreded by number of samples currentStats.normalize(currentStats.numSamples); OPENPGL_ASSERT(currentStats.isValid()); - partialWeightedMaximumAPosteriorStep(vmm, mask, currentStats, previousStats, cfg); + // Replacing the sufficient stats of the component which should not be updated with the previous suffient stats. + partialMergeSufficientStatisticsWithPriors(mask, currentStats, usePreviousStatsAsPrior, previousAsPriorMask, previousStats); + // Deriving the mixture parameter from the merged suffcient statistics using MAP and addtional global priors. + weightedMaximumAPosteriorStep(vmm, currentStats, cfg); currentEMIteration++; - // TODO: Add convergence check + + // Checking of the fitting converged (i.e., the loglikelihood only changed minimally) if (currentEMIteration > 1) { float relLogLikelihoodDifference = std::fabs(logLikelihood - previousLogLikelihood) * inv_previousLogLikelihood; @@ -925,20 +1011,196 @@ void ParallaxAwareVonMisesFisherWeightedEMFactory::partialUpda { converged = true; } - // std::cout << "logLikelihood:" << logLikelihood << "\t previousLogLikelihood: "<< previousLogLikelihood << "\t relLogLikelihoodDifference: " << - // relLogLikelihoodDifference << std::endl; previousLogLikelihood = logLikelihood; inv_previousLogLikelihood = 1.0f / std::fabs(logLikelihood); } } - previousStats += currentStats; + + // The merged sufficient stats from the last iteration are now the new previous/prior stats + previousStats = currentStats; fitStats.numSamples = numSamples; fitStats.numIterations = currentEMIteration; fitStats.summedWeightedLogLikelihood = previousLogLikelihood; - // std::cout << "converged:" << currentEMIteration << std::endl; } +#if 1 +template +void ParallaxAwareVonMisesFisherWeightedEMFactory::partialMergeSufficientStatisticsWithPriors(PartialFittingMask &mask, SufficientStatistics ¤tStats, + const bool usePreviousStatsAsPrior, + const PartialFittingMask &previousAsPriorMask, + const SufficientStatistics &previousStats) const +{ + const embree::vfloat zeros = 0.f; + + embree::vfloat sumPreviousPartialWeightsVec = 0.f; + embree::vfloat sumCurrentPartialWeightsVec = 0.f; + + float decay = 1.f; + + const int cnt = (previousStats.numComponents + VMM::VectorSize - 1) / VMM::VectorSize; + for (int k = 0; k < cnt; k++) + { + // First, check if the previous stats should be applied as prior to the current stats + if (usePreviousStatsAsPrior) + { + currentStats.sumOfWeightedDirections[k].x = + select(previousAsPriorMask.mask[k], currentStats.sumOfWeightedDirections[k].x + previousStats.sumOfWeightedDirections[k].x * decay, + currentStats.sumOfWeightedDirections[k].x); + currentStats.sumOfWeightedDirections[k].y = + select(previousAsPriorMask.mask[k], currentStats.sumOfWeightedDirections[k].y + previousStats.sumOfWeightedDirections[k].y * decay, + currentStats.sumOfWeightedDirections[k].y); + currentStats.sumOfWeightedDirections[k].z = + select(previousAsPriorMask.mask[k], currentStats.sumOfWeightedDirections[k].z + previousStats.sumOfWeightedDirections[k].z * decay, + currentStats.sumOfWeightedDirections[k].z); + + currentStats.sumOfWeightedStats[k] = + select(previousAsPriorMask.mask[k], currentStats.sumOfWeightedStats[k] + previousStats.sumOfWeightedStats[k] * decay, currentStats.sumOfWeightedStats[k]); + currentStats.sumOfDistanceWeightes[k] = select(previousAsPriorMask.mask[k], previousStats.sumOfDistanceWeightes[k] * decay, zeros); + } + // Second, reset the partial stats which sould not be updated to the prvious stats + currentStats.sumOfWeightedDirections[k].x = select(mask.mask[k], currentStats.sumOfWeightedDirections[k].x, previousStats.sumOfWeightedDirections[k].x); + currentStats.sumOfWeightedDirections[k].y = select(mask.mask[k], currentStats.sumOfWeightedDirections[k].y, previousStats.sumOfWeightedDirections[k].y); + currentStats.sumOfWeightedDirections[k].z = select(mask.mask[k], currentStats.sumOfWeightedDirections[k].z, previousStats.sumOfWeightedDirections[k].z); + + currentStats.sumOfWeightedStats[k] = select(mask.mask[k], currentStats.sumOfWeightedStats[k], previousStats.sumOfWeightedStats[k]); + currentStats.sumOfDistanceWeightes[k] = select(mask.mask[k], currentStats.sumOfDistanceWeightes[k], previousStats.sumOfDistanceWeightes[k]); + + sumPreviousPartialWeightsVec += select(mask.mask[k], previousStats.sumOfWeightedStats[k], zeros); + sumCurrentPartialWeightsVec += select(mask.mask[k], currentStats.sumOfWeightedStats[k], zeros); + } + + embree::vfloat sumTmpVec = 0.f; + // calcualting the normalization factor for the weights of the updated components + float inv_currentWeights = 1.0f; + const float sumCurrentPartialWeights = embree::reduce_add(sumCurrentPartialWeightsVec); + if (sumCurrentPartialWeights > 0.f) + { + inv_currentWeights = embree::reduce_add(sumPreviousPartialWeightsVec) / sumCurrentPartialWeights; + for (int k = 0; k < cnt; k++) + { + currentStats.sumOfWeightedStats[k] = select(mask.mask[k], currentStats.sumOfWeightedStats[k] * inv_currentWeights, currentStats.sumOfWeightedStats[k]); + + currentStats.sumOfWeightedDirections[k].x = + select(mask.mask[k], currentStats.sumOfWeightedDirections[k].x * inv_currentWeights, currentStats.sumOfWeightedDirections[k].x); + currentStats.sumOfWeightedDirections[k].y = + select(mask.mask[k], currentStats.sumOfWeightedDirections[k].y * inv_currentWeights, currentStats.sumOfWeightedDirections[k].y); + currentStats.sumOfWeightedDirections[k].z = + select(mask.mask[k], currentStats.sumOfWeightedDirections[k].z * inv_currentWeights, currentStats.sumOfWeightedDirections[k].z); + sumTmpVec += currentStats.sumOfWeightedStats[k]; + } + } + OPENPGL_ASSERT(embree::isvalid(sumTmpVec)); + float sumWTmp = embree::reduce_add(sumTmpVec); + // since only some partial components are updated + // the overall stats are the same as the previous stats + if (currentStats.normalized) + { + currentStats.sumWeights = previousStats.sumWeights; + currentStats.numSamples = sumWTmp; + } + else + { + currentStats.sumWeights = sumWTmp; + currentStats.numSamples = previousStats.numSamples; + } + + currentStats.overallNumSamples = previousStats.overallNumSamples; + + currentStats.norm = 1.f; + currentStats.inv_norm = 1.f; + if (currentStats.sumWeights > 0.f && currentStats.numSamples > 0.f) + { + currentStats.norm = currentStats.numSamples / currentStats.sumWeights; + currentStats.inv_norm = currentStats.sumWeights / currentStats.numSamples; + } + + OPENPGL_ASSERT(embree::isvalid(currentStats.norm)); + OPENPGL_ASSERT(embree::isvalid(currentStats.inv_norm)); +} +#else +template +void ParallaxAwareVonMisesFisherWeightedEMFactory::partialMergeSufficientStatisticsWithPriors(PartialFittingMask &mask, SufficientStatistics ¤tStats, + const bool usePreviousStatsAsPrior, + const PartialFittingMask &previousAsPriorMask, + const SufficientStatistics &previousStats) const +{ + const embree::vfloat zeros = 0.f; + + // the sum of the current weights of the changed components + embree::vfloat currentWeightsVec = 0.f; + // the sum of the previous weights of the unchanged components + embree::vfloat previousWeightsVec = 0.f; + // the overall sum of all previous weights + embree::vfloat sumPreviousWeightsVec = 0.f; + // embree::vfloat sumCurrentWeightsVec = 0.f; + + const int cnt = (previousStats.numComponents + VMM::VectorSize - 1) / VMM::VectorSize; + for (int k = 0; k < cnt; k++) + { + sumPreviousWeightsVec += previousStats.sumOfWeightedStats[k]; + // sumCurrentWeightsVec += currentStats.sumOfWeightedStats[k]; + + // First, check if the prior from the previous training iterations should be applied + currentStats.sumOfWeightedDirections[k].x = + select(previousAsPriorMask.mask[k], currentStats.sumOfWeightedDirections[k].x + previousStats.sumOfWeightedDirections[k].x, currentStats.sumOfWeightedDirections[k].x); + currentStats.sumOfWeightedDirections[k].y = + select(previousAsPriorMask.mask[k], currentStats.sumOfWeightedDirections[k].y + previousStats.sumOfWeightedDirections[k].y, currentStats.sumOfWeightedDirections[k].y); + currentStats.sumOfWeightedDirections[k].z = + select(previousAsPriorMask.mask[k], currentStats.sumOfWeightedDirections[k].z + previousStats.sumOfWeightedDirections[k].z, currentStats.sumOfWeightedDirections[k].z); + + currentStats.sumOfWeightedStats[k] = + select(previousAsPriorMask.mask[k], currentStats.sumOfWeightedStats[k] + previousStats.sumOfWeightedStats[k], currentStats.sumOfWeightedStats[k]); + currentStats.sumOfDistanceWeightes[k] = select(previousAsPriorMask.mask[k], previousStats.sumOfDistanceWeightes[k], zeros); + + // Second, check if the prior from the previous training iterations should be applied + currentStats.sumOfWeightedDirections[k].x = select(mask.mask[k], currentStats.sumOfWeightedDirections[k].x, previousStats.sumOfWeightedDirections[k].x); + currentStats.sumOfWeightedDirections[k].y = select(mask.mask[k], currentStats.sumOfWeightedDirections[k].y, previousStats.sumOfWeightedDirections[k].y); + currentStats.sumOfWeightedDirections[k].z = select(mask.mask[k], currentStats.sumOfWeightedDirections[k].z, previousStats.sumOfWeightedDirections[k].z); + + currentStats.sumOfWeightedStats[k] = select(mask.mask[k], currentStats.sumOfWeightedStats[k], previousStats.sumOfWeightedStats[k]); + currentStats.sumOfDistanceWeightes[k] = select(mask.mask[k], currentStats.sumOfDistanceWeightes[k], previousStats.sumOfDistanceWeightes[k]); + + currentWeightsVec += select(mask.mask[k], currentStats.sumOfWeightedStats[k], zeros); + previousWeightsVec += select(mask.mask[k], zeros, currentStats.sumOfWeightedStats[k]); + } + + float currentWeights = embree::reduce_add(currentWeightsVec); + float previousWeights = embree::reduce_add(previousWeightsVec); + float sumPreviousWeights = embree::reduce_add(sumPreviousWeightsVec); + // float sumCurrentWeights = embree::reduce_add(sumCurrentWeightsVec); + // std::cout << "sumCurrentWeights = "<< sumCurrentWeights << "\t currentStats.sumWeights = "<< currentStats.sumWeights << std::endl; + // std::cout << "sumPreviousWeights = " << sumPreviousWeights << "\t previousStats.sumWeights = "<< previousStats.sumWeights << std::endl; + + embree::vfloat sumTmpVec = 0.f; + // calcualting the normalization factor for the weights of the updated components + float inv_currentWeights = (sumPreviousWeights - previousWeights) / currentWeights; + // std::cout << "sumPreviousWeights = " << sumPreviousWeights << "\t previousWeights = "<< previousWeights << "\t currentWeights = "<< currentWeights << "\t inv_currentWeights + // = "<< inv_currentWeights << std::endl; + for (int k = 0; k < cnt; k++) + { + currentStats.sumOfWeightedStats[k] = select(mask.mask[k], currentStats.sumOfWeightedStats[k] * inv_currentWeights, currentStats.sumOfWeightedStats[k]); + currentStats.sumOfWeightedDirections[k].x = select(mask.mask[k], currentStats.sumOfWeightedDirections[k].x * inv_currentWeights, currentStats.sumOfWeightedDirections[k].x); + currentStats.sumOfWeightedDirections[k].y = select(mask.mask[k], currentStats.sumOfWeightedDirections[k].y * inv_currentWeights, currentStats.sumOfWeightedDirections[k].y); + currentStats.sumOfWeightedDirections[k].z = select(mask.mask[k], currentStats.sumOfWeightedDirections[k].z * inv_currentWeights, currentStats.sumOfWeightedDirections[k].z); + sumTmpVec += currentStats.sumOfWeightedStats[k]; + + OPENPGL_ASSERT(embree::isvalid(currentStats.sumOfWeightedStats[k])); + OPENPGL_ASSERT(embree::isvalid(currentStats.sumOfWeightedDirections[k].x)); + OPENPGL_ASSERT(embree::isvalid(currentStats.sumOfWeightedDirections[k].y)); + OPENPGL_ASSERT(embree::isvalid(currentStats.sumOfWeightedDirections[k].z)); + } + + // since only some partial components are updated + // the overall stats are the same as the previous stats + currentStats.sumWeights = previousStats.sumWeights; + currentStats.numSamples = previousStats.numSamples; + currentStats.overallNumSamples = previousStats.overallNumSamples; + currentStats.norm = previousStats.norm; + currentStats.inv_norm = previousStats.inv_norm; +} +#endif + template float ParallaxAwareVonMisesFisherWeightedEMFactory::weightedExpectationStep(VMM &vmm, SufficientStatistics &stats, UnassignedSamplesStatistics &unassignedStats, const SampleData *samples, const size_t numSamples) const @@ -962,52 +1224,52 @@ float ParallaxAwareVonMisesFisherWeightedEMFactory::weightedEx const Vector3 sampleDirection(direction.x, direction.y, direction.z); const embree::Vec3 > sampleDirectionSIMD(sampleDirection); - // check if the samples is covered by any of the components + // Calculating the soft assignment of the current sample direction for all mixture components. + // We collect the sufficient statistics for all sample directions not covered by any mixture component. if (!vmm.softAssignment(sampleDirection, softAssign)) { unassignedStats.sumOfUnassignedWeights += sampleData.weight; unassignedStats.sumUnassignedWeightedDirections += sampleDirection * sampleData.weight; -#ifdef OPENPGL_DEBUG_MODE - std::cout << "continue" << std::endl; -#endif continue; } + // Updating the sumed loglikelihood summedWeightedLogLikelihood += sampleData.weight * embree::log(softAssign.pdf); for (size_t k = 0; k < cnt; k++) { stats.sumOfWeightedDirections[k] += sampleDirectionSIMD * softAssign.assignments[k] * sampleWeight; stats.sumOfWeightedStats[k] += softAssign.assignments[k] * sampleWeight; + + OPENPGL_ASSERT(embree::isvalid(stats.sumOfWeightedDirections[k].x)); + OPENPGL_ASSERT(embree::isvalid(stats.sumOfWeightedDirections[k].y)); + OPENPGL_ASSERT(embree::isvalid(stats.sumOfWeightedDirections[k].z)); + OPENPGL_ASSERT(embree::isvalid(stats.sumOfWeightedStats[k])); } } return summedWeightedLogLikelihood; } template -void ParallaxAwareVonMisesFisherWeightedEMFactory::estimateMAPWeights(VMM &vmm, const SufficientStatistics ¤tStats, - const SufficientStatistics &previousStats, const float &_weightPrior) const +void ParallaxAwareVonMisesFisherWeightedEMFactory::estimateMAPWeights(VMM &vmm, const SufficientStatistics ¤tStats, const float &_weightPrior) const { const int cnt = (vmm._numComponents + VMM::VectorSize - 1) / VMM::VectorSize; const size_t numComponents = vmm._numComponents; - const embree::vfloat weightPrior(_weightPrior); - - const embree::vfloat numSamples = currentStats.numSamples + previousStats.numSamples; - // const vfloat numSamples = currentStats.numSamples + previousStats.overallNumSamples; + const embree::vfloat numSamples = currentStats.numSamples; for (size_t k = 0; k < cnt; k++) { - //_sumWeights += currentStats.sumOfWeightedStats[k]; - embree::vfloat weight = (currentStats.sumOfWeightedStats[k] + previousStats.sumOfWeightedStats[k]); + OPENPGL_ASSERT(embree::isvalid(vmm._weights[k])); + embree::vfloat weight = (currentStats.sumOfWeightedStats[k]); weight = (weightPrior + (weight)) / ((weightPrior * numComponents) + numSamples); - // vfloat weight = ( currentStats.sumOfWeightedStats[k]/* + previousStats.sumOfWeightedStats[k]*/ ) / ( sumWeights ); - // weight = ( weightPrior + ( weight * numSamples ) ) / (( weightPrior * numComponents ) + numSamples ); vmm._weights[k] = weight; + OPENPGL_ASSERT(embree::isvalid(vmm._weights[k])); } // TODO: find better more efficient way + // Ensuring that the weights for unused SIMD vector entries are zero if (vmm._numComponents % VMM::VectorSize > 0) { for (size_t i = vmm._numComponents % VMM::VectorSize; i < VMM::VectorSize; i++) @@ -1019,16 +1281,9 @@ void ParallaxAwareVonMisesFisherWeightedEMFactory::estimateMAP template void ParallaxAwareVonMisesFisherWeightedEMFactory::estimateMAPMeanDirectionAndConcentration(VMM &vmm, const SufficientStatistics ¤tStats, - const SufficientStatistics &previousStats, const Configuration &cfg) const { - const embree::vfloat currentNumSamples = currentStats.numSamples; - const embree::vfloat previousNumSamples = previousStats.numSamples; - const embree::vfloat numSamples = currentNumSamples + previousNumSamples; - const embree::vfloat overallNumSamples = currentStats.numSamples + previousStats.overallNumSamples; - - const embree::vfloat currentEstimationWeight = currentNumSamples / numSamples; - const embree::vfloat previousEstimationWeight = 1.0f - currentEstimationWeight; + const embree::vfloat overallNumSamples = currentStats.overallNumSamples; const embree::vfloat meanCosinePrior = cfg.meanCosinePrior; const embree::vfloat meanCosinePriorStrength = cfg.meanCosinePriorStrength; @@ -1038,20 +1293,11 @@ void ParallaxAwareVonMisesFisherWeightedEMFactory::estimateMAP for (size_t k = 0; k < cnt; k++) { - // const vfloat partialNumSamples = vmm._weights[k] * numSamples; const embree::vfloat partialNumSamples = vmm._weights[k] * overallNumSamples; - embree::Vec3 > currentMeanDirection; - currentMeanDirection.x = select(currentStats.sumOfWeightedStats[k] > 0.0f, currentStats.sumOfWeightedDirections[k].x / currentStats.sumOfWeightedStats[k], 0.0f); - currentMeanDirection.y = select(currentStats.sumOfWeightedStats[k] > 0.0f, currentStats.sumOfWeightedDirections[k].y / currentStats.sumOfWeightedStats[k], 0.0f); - currentMeanDirection.z = select(currentStats.sumOfWeightedStats[k] > 0.0f, currentStats.sumOfWeightedDirections[k].z / currentStats.sumOfWeightedStats[k], 0.0f); - - // TODO: find a better design to precompute the previousMeanDirection - embree::Vec3 > previousMeanDirection; - previousMeanDirection.x = select(previousStats.sumOfWeightedStats[k] > 0.0f, previousStats.sumOfWeightedDirections[k].x / previousStats.sumOfWeightedStats[k], 0.0f); - previousMeanDirection.y = select(previousStats.sumOfWeightedStats[k] > 0.0f, previousStats.sumOfWeightedDirections[k].y / previousStats.sumOfWeightedStats[k], 0.0f); - previousMeanDirection.z = select(previousStats.sumOfWeightedStats[k] > 0.0f, previousStats.sumOfWeightedDirections[k].z / previousStats.sumOfWeightedStats[k], 0.0f); - - embree::Vec3 > meanDirection = currentMeanDirection * currentEstimationWeight + previousMeanDirection * previousEstimationWeight; + embree::Vec3 > meanDirection; + meanDirection.x = select(currentStats.sumOfWeightedStats[k] > 0.0f, currentStats.sumOfWeightedDirections[k].x / currentStats.sumOfWeightedStats[k], 0.0f); + meanDirection.y = select(currentStats.sumOfWeightedStats[k] > 0.0f, currentStats.sumOfWeightedDirections[k].y / currentStats.sumOfWeightedStats[k], 0.0f); + meanDirection.z = select(currentStats.sumOfWeightedStats[k] > 0.0f, currentStats.sumOfWeightedDirections[k].z / currentStats.sumOfWeightedStats[k], 0.0f); embree::vfloat meanCosine = length(meanDirection); @@ -1059,14 +1305,19 @@ void ParallaxAwareVonMisesFisherWeightedEMFactory::estimateMAP vmm._meanDirections[k].y = select(meanCosine > 0.0f, meanDirection.y / meanCosine, vmm._meanDirections[k].y); vmm._meanDirections[k].z = select(meanCosine > 0.0f, meanDirection.z / meanCosine, vmm._meanDirections[k].z); + // Applying a global prior to the mean cosines of the mixture components. + // The strength of the prior depends on the partial number of samples associated with each component. meanCosine = (meanCosinePrior * meanCosinePriorStrength + meanCosine * partialNumSamples) / (meanCosinePriorStrength + partialNumSamples); meanCosine = embree::min(maxMeanCosine, meanCosine); vmm._meanCosines[k] = meanCosine; vmm._kappas[k] = MeanCosineToKappa >(meanCosine); + OPENPGL_ASSERT(embree::isvalid(vmm._meanCosines[k])); + OPENPGL_ASSERT(embree::isvalid(vmm._kappas[k])); } // TODO: find better more efficient way + // Ensuring that the component parameters for unused SIMD vector entries are cleared if (rem > 0) { for (size_t i = rem; i < VMM::VectorSize; i++) @@ -1088,140 +1339,13 @@ void ParallaxAwareVonMisesFisherWeightedEMFactory::estimateMAP template void ParallaxAwareVonMisesFisherWeightedEMFactory::weightedMaximumAPosteriorStep(VMM &vmm, const SufficientStatistics ¤tStats, - const SufficientStatistics &previousStats, const Configuration &cfg) const -{ - // Estimating components weights - estimateMAPWeights(vmm, currentStats, previousStats, cfg.weightPrior); - - // Estimating mean and concentration - estimateMAPMeanDirectionAndConcentration(vmm, currentStats, previousStats, cfg); -} - -template -void ParallaxAwareVonMisesFisherWeightedEMFactory::partialWeightedMaximumAPosteriorStep(VMM &vmm, const PartialFittingMask &mask, - SufficientStatistics ¤tStats, SufficientStatistics &previousStats, - const Configuration &cfg) const + const Configuration &cfg) const { - // Estimating components weights - estimatePartialMAPWeights(vmm, mask, currentStats, previousStats, cfg.weightPrior); - - // Estimating mean and concentration - estimatePartialMAPMeanDirectionAndConcentration(vmm, mask, currentStats, previousStats, cfg); -} - -template -void ParallaxAwareVonMisesFisherWeightedEMFactory::estimatePartialMAPWeights(VMM &vmm, const PartialFittingMask &mask, SufficientStatistics ¤tStats, - SufficientStatistics &previousStats, const float &_weightPrior) const -{ - const embree::vfloat zeros(0.0f); - const int cnt = (vmm._numComponents + VMM::VectorSize - 1) / VMM::VectorSize; - - const size_t numComponents = vmm._numComponents; - - const embree::vfloat weightPrior(_weightPrior); - - const embree::vfloat numSamples = currentStats.numSamples + previousStats.numSamples; - // const vfloat numSamples = currentStats.numSamples + previousStats.overallNumSamples; - - embree::vfloat sumWeights(0.0f); - embree::vfloat sumPartialWeights(0.0f); + // MAP estimation of the components weights + estimateMAPWeights(vmm, currentStats, cfg.weightPrior); - for (size_t k = 0; k < cnt; k++) - { - //_sumWeights += currentStats.sumOfWeightedStats[k]; - embree::vfloat weight = (currentStats.sumOfWeightedStats[k] + previousStats.sumOfWeightedStats[k]); - weight = (weightPrior + (weight)) / ((weightPrior * numComponents) + numSamples); - // vfloat weight = ( currentStats.sumOfWeightedStats[k]/* + previousStats.sumOfWeightedStats[k]*/ ) / ( sumWeights ); - // weight = ( weightPrior + ( weight * numSamples ) ) / (( weightPrior * numComponents ) + numSamples ); - - sumPartialWeights += select(mask.mask[k], weight, zeros); - sumWeights += select(mask.mask[k], zeros, vmm._weights[k]); - - vmm._weights[k] = select(mask.mask[k], weight, vmm._weights[k]); - } - - embree::vfloat inv_sumPartialWeights = 1.0f / reduce_add(sumPartialWeights); - inv_sumPartialWeights *= 1.0f - reduce_add(sumWeights); - for (size_t k = 0; k < cnt; k++) - { - vmm._weights[k] = select(mask.mask[k], vmm._weights[k] * inv_sumPartialWeights, vmm._weights[k]); - } - - // TODO: find better more efficient way - if (vmm._numComponents % VMM::VectorSize > 0) - { - for (size_t i = vmm._numComponents % VMM::VectorSize; i < VMM::VectorSize; i++) - { - vmm._weights[cnt - 1][i] = 0.0f; - } - } -} - -template -void ParallaxAwareVonMisesFisherWeightedEMFactory::estimatePartialMAPMeanDirectionAndConcentration(VMM &vmm, const PartialFittingMask &mask, - SufficientStatistics ¤tStats, - SufficientStatistics &previousStats, - const Configuration &cfg) const -{ - const embree::vfloat currentNumSamples = currentStats.numSamples; - const embree::vfloat previousNumSamples = previousStats.numSamples; - const embree::vfloat numSamples = currentNumSamples + previousNumSamples; - const embree::vfloat overallNumSamples = currentStats.numSamples + previousStats.overallNumSamples; - - const embree::vfloat currentEstimationWeight = currentNumSamples / numSamples; - const embree::vfloat previousEstimationWeight = 1.0f - currentEstimationWeight; - - const embree::vfloat meanCosinePrior = cfg.meanCosinePrior; - const embree::vfloat meanCosinePriorStrength = cfg.meanCosinePriorStrength; - const embree::vfloat maxMeanCosine = cfg.maxMeanCosine; - const int cnt = (vmm._numComponents + VMM::VectorSize - 1) / VMM::VectorSize; - const int rem = vmm._numComponents % VMM::VectorSize; - - for (size_t k = 0; k < cnt; k++) - { - // const vfloat partialNumSamples = vmm._weights[k] * numSamples; - const embree::vfloat partialNumSamples = vmm._weights[k] * overallNumSamples; - embree::Vec3 > currentMeanDirection; - currentMeanDirection.x = select(currentStats.sumOfWeightedStats[k] > 0.0f, currentStats.sumOfWeightedDirections[k].x / currentStats.sumOfWeightedStats[k], 0.0f); - currentMeanDirection.y = select(currentStats.sumOfWeightedStats[k] > 0.0f, currentStats.sumOfWeightedDirections[k].y / currentStats.sumOfWeightedStats[k], 0.0f); - currentMeanDirection.z = select(currentStats.sumOfWeightedStats[k] > 0.0f, currentStats.sumOfWeightedDirections[k].z / currentStats.sumOfWeightedStats[k], 0.0f); - - // TODO: find a better design to precompute the previousMeanDirection - embree::Vec3 > previousMeanDirection; - previousMeanDirection.x = select(previousStats.sumOfWeightedStats[k] > 0.0f, previousStats.sumOfWeightedDirections[k].x / previousStats.sumOfWeightedStats[k], 0.0f); - previousMeanDirection.y = select(previousStats.sumOfWeightedStats[k] > 0.0f, previousStats.sumOfWeightedDirections[k].y / previousStats.sumOfWeightedStats[k], 0.0f); - previousMeanDirection.z = select(previousStats.sumOfWeightedStats[k] > 0.0f, previousStats.sumOfWeightedDirections[k].z / previousStats.sumOfWeightedStats[k], 0.0f); - - embree::Vec3 > meanDirection = currentMeanDirection * currentEstimationWeight + previousMeanDirection * previousEstimationWeight; - - embree::vfloat meanCosine = embree::length(meanDirection); - - vmm._meanDirections[k].x = select(mask.mask[k], select(meanCosine > 0.0f, meanDirection.x / meanCosine, vmm._meanDirections[k].x), vmm._meanDirections[k].x); - vmm._meanDirections[k].y = select(mask.mask[k], select(meanCosine > 0.0f, meanDirection.y / meanCosine, vmm._meanDirections[k].y), vmm._meanDirections[k].y); - vmm._meanDirections[k].z = select(mask.mask[k], select(meanCosine > 0.0f, meanDirection.z / meanCosine, vmm._meanDirections[k].z), vmm._meanDirections[k].z); - - meanCosine = (meanCosinePrior * meanCosinePriorStrength + meanCosine * partialNumSamples) / (meanCosinePriorStrength + partialNumSamples); - - meanCosine = embree::min(maxMeanCosine, meanCosine); - vmm._meanCosines[k] = select(mask.mask[k], meanCosine, vmm._meanCosines[k]); - vmm._kappas[k] = select(mask.mask[k], MeanCosineToKappa >(meanCosine), vmm._kappas[k]); - } - - // TODO: find better more efficient way - if (rem > 0) - { - for (size_t i = rem; i < VMM::VectorSize; i++) - { - vmm._meanDirections[cnt - 1].x[i] = 0.0f; - vmm._meanDirections[cnt - 1].y[i] = 0.0f; - vmm._meanDirections[cnt - 1].z[i] = 1.0f; - - vmm._meanCosines[cnt - 1][i] = 0.0f; - vmm._kappas[cnt - 1][i] = 0.0f; - } - } - - vmm._calculateNormalization(); + // MAP estimation of the mean and concentration + estimateMAPMeanDirectionAndConcentration(vmm, currentStats, cfg); } template @@ -1666,6 +1790,16 @@ void ParallaxAwareVonMisesFisherWeightedEMFactory::PartialFitt } } +template +void ParallaxAwareVonMisesFisherWeightedEMFactory::PartialFittingMask::resetToTrue() +{ + const embree::vbool vTrue(true); + for (size_t k = 0; k < ((VMM::MaxComponents + (VMM::VectorSize - 1)) / VMM::VectorSize); k++) + { + mask[k] = vTrue; + } +} + template void ParallaxAwareVonMisesFisherWeightedEMFactory::PartialFittingMask::resetToTrue(const size_t &numComponents) { diff --git a/openpgl/directional/vmm/VMMChiSquareComponentMerger.h b/openpgl/directional/vmm/VMMChiSquareComponentMerger.h index d27e6787..4bbaa6a3 100644 --- a/openpgl/directional/vmm/VMMChiSquareComponentMerger.h +++ b/openpgl/directional/vmm/VMMChiSquareComponentMerger.h @@ -10,26 +10,28 @@ namespace openpgl { -template +template struct VonMisesFisherChiSquareComponentMerger { public: typedef typename TVMMFactory::Distribution VMM; typedef typename TVMMFactory::SufficientStatistics SufficientStatistics; - typedef typename VonMisesFisherChiSquareComponentSplitter::ComponentSplitStatistics ComponentSplitStatistics; + typedef typename TSplitter::ComponentSplitStatistics ComponentSplitStatistics; // typedef std::integral_constant NumVectors; float MergeNext(VMM &vmm) const; bool ThresholdedMergeNext(VMM &vmm, const float &mergeThreshold, float &mergeCost) const; - bool ThresholdedMergeNext(VMM &vmm, const float &mergeThreshold, float &mergeCost, SufficientStatistics &suffStats, ComponentSplitStatistics &splitStats) const; + bool ThresholdedMergeNext(VMM &vmm, const float mergeThreshold, const float splitThreshold, const bool considerSplitThreshold, float &mergeCost, + SufficientStatistics &suffStats, ComponentSplitStatistics &splitStats) const; float CalculateMergeCost(const VMM &vmm, const size_t &idx0, const size_t &idx1) const; size_t PerformMerging(VMM &vmm, const float &mergeThreshold) const; - size_t PerformMerging(VMM &vmm, const float &mergeThreshold, SufficientStatistics &suffStats, ComponentSplitStatistics &splitStats) const; + size_t PerformMerging(VMM &vmm, const float mergeThreshold, const float splitThreshold, const bool considerSplitThreshold, SufficientStatistics &suffStats, + ComponentSplitStatistics &splitStats) const; private: inline float _IntegratedProduct(const Vector3 &meanDirection0, const float &kappa0, const float &normalization0, const Vector3 &meanDirection1, const float &kappa1, @@ -42,8 +44,8 @@ struct VonMisesFisherChiSquareComponentMerger const float &normalization1, Vector3 &meanDirection, float &kappa, float &normalization) const; }; -template -size_t VonMisesFisherChiSquareComponentMerger::PerformMerging(VMM &vmm, const float &mergeThreshold) const +template +size_t VonMisesFisherChiSquareComponentMerger::PerformMerging(VMM &vmm, const float &mergeThreshold) const { bool stopMerging = false; size_t totalMergeCount = 0; @@ -63,9 +65,10 @@ size_t VonMisesFisherChiSquareComponentMerger::PerformMerging(VMM & return totalMergeCount; } -template -size_t VonMisesFisherChiSquareComponentMerger::PerformMerging(VMM &vmm, const float &mergeThreshold, SufficientStatistics &suffStats, - ComponentSplitStatistics &splitStats) const +template +size_t VonMisesFisherChiSquareComponentMerger::PerformMerging(VMM &vmm, const float mergeThreshold, const float splitThreshold, + const bool considerSplitThreshold, SufficientStatistics &suffStats, + ComponentSplitStatistics &splitStats) const { bool stopMerging = false; size_t totalMergeCount = 0; @@ -74,7 +77,7 @@ size_t VonMisesFisherChiSquareComponentMerger::PerformMerging(VMM & { float mergeCost = 0.0f; stopMerging = true; - bool mergedComponents = ThresholdedMergeNext(vmm, mergeThreshold, mergeCost, suffStats, splitStats); + bool mergedComponents = ThresholdedMergeNext(vmm, mergeThreshold, splitThreshold, considerSplitThreshold, mergeCost, suffStats, splitStats); stopMerging = !mergedComponents; // std::cout << "merge: " << "\tmergedComponents: " << mergedComponents << "\tmergeCost: " << mergeCost << std::endl; @@ -89,8 +92,8 @@ size_t VonMisesFisherChiSquareComponentMerger::PerformMerging(VMM & return totalMergeCount; } -template -float VonMisesFisherChiSquareComponentMerger::CalculateMergeCost(const VMM &vmm, const size_t &idx0, const size_t &idx1) const +template +float VonMisesFisherChiSquareComponentMerger::CalculateMergeCost(const VMM &vmm, const size_t &idx0, const size_t &idx1) const { const div_t div0 = div(idx0, VMM::VectorSize); float weight0 = vmm._weights[div0.quot][div0.rem]; @@ -165,8 +168,8 @@ float VonMisesFisherChiSquareComponentMerger::CalculateMergeCost(co return chiSquareIJ; } -template -bool VonMisesFisherChiSquareComponentMerger::ThresholdedMergeNext(VMM &vmm, const float &mergeThreshold, float &mergeCost) const +template +bool VonMisesFisherChiSquareComponentMerger::ThresholdedMergeNext(VMM &vmm, const float &mergeThreshold, float &mergeCost) const { // std::cout << vmm.toString()<::ThresholdedMergeNext(V return foundMergeCandidates; } -template -bool VonMisesFisherChiSquareComponentMerger::ThresholdedMergeNext(VMM &vmm, const float &mergeThreshold, float &mergeCost, SufficientStatistics &suffStats, - ComponentSplitStatistics &splitStats) const +template +bool VonMisesFisherChiSquareComponentMerger::ThresholdedMergeNext(VMM &vmm, const float mergeThreshold, const float splitThreshold, + const bool considerSplitThreshold, float &mergeCost, SufficientStatistics &suffStats, + ComponentSplitStatistics &splitStats) const { OPENPGL_ASSERT(splitStats.isValid()); int K = vmm._numComponents; @@ -219,7 +223,9 @@ bool VonMisesFisherChiSquareComponentMerger::ThresholdedMergeNext(V { const div_t tmpJ = div(j, static_cast(VMM::VectorSize)); float mergeCost = CalculateMergeCost(vmm, i, j); - if (mergeCost < mergeThreshold && mergeCost < minMergeCost && splitStats.numSamples[tmpI.quot][tmpI.rem] > 0.0f && splitStats.numSamples[tmpJ.quot][tmpJ.rem] > 0.0f) + if (mergeCost < mergeThreshold && mergeCost < minMergeCost && splitStats.numSamples[tmpI.quot][tmpI.rem] > 0.0f && splitStats.numSamples[tmpJ.quot][tmpJ.rem] > 0.0f && + (!considerSplitThreshold || + (considerSplitThreshold && (splitStats.chiSquareMCEstimates[tmpI.quot][tmpI.rem] + splitStats.chiSquareMCEstimates[tmpJ.quot][tmpJ.rem] < splitThreshold)))) { mergeCandidateI = i; mergeCandidateJ = j; @@ -262,8 +268,8 @@ bool VonMisesFisherChiSquareComponentMerger::ThresholdedMergeNext(V return foundMergeCandidates; } -template -float VonMisesFisherChiSquareComponentMerger::MergeNext(VMM &vmm) const +template +float VonMisesFisherChiSquareComponentMerger::MergeNext(VMM &vmm) const { int K = vmm._numComponents; int mergeCandidateI = 0; @@ -288,9 +294,9 @@ float VonMisesFisherChiSquareComponentMerger::MergeNext(VMM &vmm) c return minMergeCost; } -template -float VonMisesFisherChiSquareComponentMerger::_IntegratedProduct(const Vector3 &meanDirection0, const float &kappa0, const float &normalization0, - const Vector3 &meanDirection1, const float &kappa1, const float &normalization1) const +template +float VonMisesFisherChiSquareComponentMerger::_IntegratedProduct(const Vector3 &meanDirection0, const float &kappa0, const float &normalization0, + const Vector3 &meanDirection1, const float &kappa1, const float &normalization1) const { Vector3 productMeanDirection = kappa0 * meanDirection0 + kappa1 * meanDirection1; @@ -318,10 +324,10 @@ float VonMisesFisherChiSquareComponentMerger::_IntegratedProduct(co return scale; } -template -float VonMisesFisherChiSquareComponentMerger::_IntegratedDivision(const Vector3 &meanDirection0, const float &kappa0, const float &normalization0, - const Vector3 &meanDirection1, const float &kappa1, const float &normalization1, - const float &eMinus2Kappa1) const +template +float VonMisesFisherChiSquareComponentMerger::_IntegratedDivision(const Vector3 &meanDirection0, const float &kappa0, const float &normalization0, + const Vector3 &meanDirection1, const float &kappa1, const float &normalization1, + const float &eMinus2Kappa1) const { Vector3 productMeanDirection = kappa0 * meanDirection0 + kappa1 * meanDirection1; @@ -350,10 +356,10 @@ float VonMisesFisherChiSquareComponentMerger::_IntegratedDivision(c return scale; } -template -float VonMisesFisherChiSquareComponentMerger::_Product(const Vector3 &meanDirection0, const float &kappa0, const float &normalization0, const Vector3 &meanDirection1, - const float &kappa1, const float &normalization1, Vector3 &productMeanDirection, float &productKappa, - float &productNormalization) const +template +float VonMisesFisherChiSquareComponentMerger::_Product(const Vector3 &meanDirection0, const float &kappa0, const float &normalization0, + const Vector3 &meanDirection1, const float &kappa1, const float &normalization1, + Vector3 &productMeanDirection, float &productKappa, float &productNormalization) const { productMeanDirection = kappa0 * meanDirection0 + kappa1 * meanDirection1; productKappa = embree::sqrt(dot(productMeanDirection, productMeanDirection)); diff --git a/openpgl/directional/vmm/VMMChiSquareComponentSplitter.h b/openpgl/directional/vmm/VMMChiSquareComponentSplitter.h index ef12ddb8..9f4f0d1a 100644 --- a/openpgl/directional/vmm/VMMChiSquareComponentSplitter.h +++ b/openpgl/directional/vmm/VMMChiSquareComponentSplitter.h @@ -278,7 +278,7 @@ void VonMisesFisherChiSquareComponentSplitter::PerformSplitting(VMM // std::cout << "suffStatistics: " << suffStatistics.toString() << std::endl; if (doPartialRefit) { - vmmFactory.partialUpdateMixture(vmm, mask, suffStatistics, data, numData, factoryCfg, vmmFitStats); + vmmFactory.partialUpdateMixture(vmm, mask, suffStatistics, true, data, numData, factoryCfg, vmmFitStats); // std::cout << "vmmpartialUpdate: " << vmm.toString() << std::endl; splitItr++; } @@ -300,8 +300,9 @@ void VonMisesFisherChiSquareComponentSplitter::PerformRecursiveSpli const size_t &numData, const typename VMMFactory::Configuration factoryCfg) const { PartialFittingMask mask; + PartialFittingMask previousAsPriorMask; + previousAsPriorMask.resetToFalse(); ComponentSplitStatistics splitStatistics; - SufficientStatistics tempSuffStatistics = suffStatistics; // bool stopSplitting = false; // size_t splitItr = 0; @@ -330,14 +331,14 @@ void VonMisesFisherChiSquareComponentSplitter::PerformRecursiveSpli if (splitComps[k].chiSquareEst > splitThreshold && vmm._numComponents < VMM::MaxComponents) { #ifndef OPENPGL_USE_THREE_SPLIT - bool splitSuccess = SplitComponent(vmm, splitStatistics, tempSuffStatistics, splitComps[k].componentIndex); + bool splitSuccess = SplitComponent(vmm, splitStatistics, suffStatistics, splitComps[k].componentIndex); if (splitSuccess) { mask.setToTrue(splitComps[k].componentIndex); mask.setToTrue(vmm._numComponents - 1); } #else - bool splitSuccess = SplitComponentIntoThree(vmm, splitStatistics, tempSuffStatistics, splitComps[k].componentIndex); + bool splitSuccess = SplitComponentIntoThree(vmm, splitStatistics, suffStatistics, splitComps[k].componentIndex); if (splitSuccess) { mask.setToTrue(splitComps[k].componentIndex); @@ -357,11 +358,7 @@ void VonMisesFisherChiSquareComponentSplitter::PerformRecursiveSpli } if (numSplits > 0) { - tempSuffStatistics.clear(vmm._numComponents); - vmmFactory.partialUpdateMixture(vmm, mask, tempSuffStatistics, data, numData, factoryCfg, vmmFitStats); - // std::cout << "tempSuffStatistics" << std::endl << tempSuffStatistics.toString() << std::endl; - suffStatistics.setNumComponents(vmm._numComponents); - suffStatistics.maskedReplace(mask, tempSuffStatistics); + vmmFactory.partialUpdateMixture(vmm, mask, false, previousAsPriorMask, suffStatistics, data, numData, factoryCfg, vmmFitStats); } // std::cout << "vmmpartialUpdate: " << vmm.toString() << std::endl; // splitItr++; From 2ffaa3db5ebe440a0f633d8d882bd4645dd73c47 Mon Sep 17 00:00:00 2001 From: Sebastian Date: Tue, 20 Oct 2026 09:16:45 +0200 Subject: [PATCH 5/5] [sam] Some additional documentation and cleanups --- ...llaxAwareVonMisesFisherWeightedEMFactory.h | 101 +++--------------- 1 file changed, 13 insertions(+), 88 deletions(-) diff --git a/openpgl/directional/vmm/ParallaxAwareVonMisesFisherWeightedEMFactory.h b/openpgl/directional/vmm/ParallaxAwareVonMisesFisherWeightedEMFactory.h index 1e85d95d..4cf7c587 100644 --- a/openpgl/directional/vmm/ParallaxAwareVonMisesFisherWeightedEMFactory.h +++ b/openpgl/directional/vmm/ParallaxAwareVonMisesFisherWeightedEMFactory.h @@ -726,7 +726,6 @@ void ParallaxAwareVonMisesFisherWeightedEMFactory::SufficientS sumOfWeightedStats[k] *= normVec; } normalized = true; - // TODO: check this seems to be wrong overallNumSamples = numSamples; } /* */ @@ -1024,7 +1023,6 @@ void ParallaxAwareVonMisesFisherWeightedEMFactory::partialUpda fitStats.summedWeightedLogLikelihood = previousLogLikelihood; } -#if 1 template void ParallaxAwareVonMisesFisherWeightedEMFactory::partialMergeSufficientStatisticsWithPriors(PartialFittingMask &mask, SufficientStatistics ¤tStats, const bool usePreviousStatsAsPrior, @@ -1033,6 +1031,8 @@ void ParallaxAwareVonMisesFisherWeightedEMFactory::partialMerg { const embree::vfloat zeros = 0.f; + // The old (previous) and new (current) sum of weights for the partial updated components. + // In the end the weights of the partial updated components must be similar embree::vfloat sumPreviousPartialWeightsVec = 0.f; embree::vfloat sumCurrentPartialWeightsVec = 0.f; @@ -1044,6 +1044,7 @@ void ParallaxAwareVonMisesFisherWeightedEMFactory::partialMerg // First, check if the previous stats should be applied as prior to the current stats if (usePreviousStatsAsPrior) { + // Merge the current and the previous stats if allowed by the previousAsPriorMask. currentStats.sumOfWeightedDirections[k].x = select(previousAsPriorMask.mask[k], currentStats.sumOfWeightedDirections[k].x + previousStats.sumOfWeightedDirections[k].x * decay, currentStats.sumOfWeightedDirections[k].x); @@ -1058,24 +1059,28 @@ void ParallaxAwareVonMisesFisherWeightedEMFactory::partialMerg select(previousAsPriorMask.mask[k], currentStats.sumOfWeightedStats[k] + previousStats.sumOfWeightedStats[k] * decay, currentStats.sumOfWeightedStats[k]); currentStats.sumOfDistanceWeightes[k] = select(previousAsPriorMask.mask[k], previousStats.sumOfDistanceWeightes[k] * decay, zeros); } - // Second, reset the partial stats which sould not be updated to the prvious stats + // Second, reset the stats that should not be updated by the partial fit. + // These stats are reseted to the ones from the previous stats. currentStats.sumOfWeightedDirections[k].x = select(mask.mask[k], currentStats.sumOfWeightedDirections[k].x, previousStats.sumOfWeightedDirections[k].x); currentStats.sumOfWeightedDirections[k].y = select(mask.mask[k], currentStats.sumOfWeightedDirections[k].y, previousStats.sumOfWeightedDirections[k].y); currentStats.sumOfWeightedDirections[k].z = select(mask.mask[k], currentStats.sumOfWeightedDirections[k].z, previousStats.sumOfWeightedDirections[k].z); + // Update the sum of weights and wieghted distances for the partial updated components. currentStats.sumOfWeightedStats[k] = select(mask.mask[k], currentStats.sumOfWeightedStats[k], previousStats.sumOfWeightedStats[k]); currentStats.sumOfDistanceWeightes[k] = select(mask.mask[k], currentStats.sumOfDistanceWeightes[k], previousStats.sumOfDistanceWeightes[k]); + // Calculate the previous and current sum of weights and distances for the partial updated components. sumPreviousPartialWeightsVec += select(mask.mask[k], previousStats.sumOfWeightedStats[k], zeros); sumCurrentPartialWeightsVec += select(mask.mask[k], currentStats.sumOfWeightedStats[k], zeros); } embree::vfloat sumTmpVec = 0.f; - // calcualting the normalization factor for the weights of the updated components + // Calcualting the normalization factor for the weights of the updated components float inv_currentWeights = 1.0f; const float sumCurrentPartialWeights = embree::reduce_add(sumCurrentPartialWeightsVec); if (sumCurrentPartialWeights > 0.f) { + // Correction factor so that the sum of weights for the updated components stays the same as before. inv_currentWeights = embree::reduce_add(sumPreviousPartialWeightsVec) / sumCurrentPartialWeights; for (int k = 0; k < cnt; k++) { @@ -1092,8 +1097,8 @@ void ParallaxAwareVonMisesFisherWeightedEMFactory::partialMerg } OPENPGL_ASSERT(embree::isvalid(sumTmpVec)); float sumWTmp = embree::reduce_add(sumTmpVec); - // since only some partial components are updated - // the overall stats are the same as the previous stats + // Since only some partial components are updated the overall stats for numSamples, sumWeights, ..., + // should not change. Nevertheless, due to some numerical differences we recalculate the stats. if (currentStats.normalized) { currentStats.sumWeights = previousStats.sumWeights; @@ -1105,8 +1110,10 @@ void ParallaxAwareVonMisesFisherWeightedEMFactory::partialMerg currentStats.numSamples = previousStats.numSamples; } + // The overall number of samples stays the same. currentStats.overallNumSamples = previousStats.overallNumSamples; + // Updating/Recalculating the weights to num samples normalization factor and its inverse. currentStats.norm = 1.f; currentStats.inv_norm = 1.f; if (currentStats.sumWeights > 0.f && currentStats.numSamples > 0.f) @@ -1118,88 +1125,6 @@ void ParallaxAwareVonMisesFisherWeightedEMFactory::partialMerg OPENPGL_ASSERT(embree::isvalid(currentStats.norm)); OPENPGL_ASSERT(embree::isvalid(currentStats.inv_norm)); } -#else -template -void ParallaxAwareVonMisesFisherWeightedEMFactory::partialMergeSufficientStatisticsWithPriors(PartialFittingMask &mask, SufficientStatistics ¤tStats, - const bool usePreviousStatsAsPrior, - const PartialFittingMask &previousAsPriorMask, - const SufficientStatistics &previousStats) const -{ - const embree::vfloat zeros = 0.f; - - // the sum of the current weights of the changed components - embree::vfloat currentWeightsVec = 0.f; - // the sum of the previous weights of the unchanged components - embree::vfloat previousWeightsVec = 0.f; - // the overall sum of all previous weights - embree::vfloat sumPreviousWeightsVec = 0.f; - // embree::vfloat sumCurrentWeightsVec = 0.f; - - const int cnt = (previousStats.numComponents + VMM::VectorSize - 1) / VMM::VectorSize; - for (int k = 0; k < cnt; k++) - { - sumPreviousWeightsVec += previousStats.sumOfWeightedStats[k]; - // sumCurrentWeightsVec += currentStats.sumOfWeightedStats[k]; - - // First, check if the prior from the previous training iterations should be applied - currentStats.sumOfWeightedDirections[k].x = - select(previousAsPriorMask.mask[k], currentStats.sumOfWeightedDirections[k].x + previousStats.sumOfWeightedDirections[k].x, currentStats.sumOfWeightedDirections[k].x); - currentStats.sumOfWeightedDirections[k].y = - select(previousAsPriorMask.mask[k], currentStats.sumOfWeightedDirections[k].y + previousStats.sumOfWeightedDirections[k].y, currentStats.sumOfWeightedDirections[k].y); - currentStats.sumOfWeightedDirections[k].z = - select(previousAsPriorMask.mask[k], currentStats.sumOfWeightedDirections[k].z + previousStats.sumOfWeightedDirections[k].z, currentStats.sumOfWeightedDirections[k].z); - - currentStats.sumOfWeightedStats[k] = - select(previousAsPriorMask.mask[k], currentStats.sumOfWeightedStats[k] + previousStats.sumOfWeightedStats[k], currentStats.sumOfWeightedStats[k]); - currentStats.sumOfDistanceWeightes[k] = select(previousAsPriorMask.mask[k], previousStats.sumOfDistanceWeightes[k], zeros); - - // Second, check if the prior from the previous training iterations should be applied - currentStats.sumOfWeightedDirections[k].x = select(mask.mask[k], currentStats.sumOfWeightedDirections[k].x, previousStats.sumOfWeightedDirections[k].x); - currentStats.sumOfWeightedDirections[k].y = select(mask.mask[k], currentStats.sumOfWeightedDirections[k].y, previousStats.sumOfWeightedDirections[k].y); - currentStats.sumOfWeightedDirections[k].z = select(mask.mask[k], currentStats.sumOfWeightedDirections[k].z, previousStats.sumOfWeightedDirections[k].z); - - currentStats.sumOfWeightedStats[k] = select(mask.mask[k], currentStats.sumOfWeightedStats[k], previousStats.sumOfWeightedStats[k]); - currentStats.sumOfDistanceWeightes[k] = select(mask.mask[k], currentStats.sumOfDistanceWeightes[k], previousStats.sumOfDistanceWeightes[k]); - - currentWeightsVec += select(mask.mask[k], currentStats.sumOfWeightedStats[k], zeros); - previousWeightsVec += select(mask.mask[k], zeros, currentStats.sumOfWeightedStats[k]); - } - - float currentWeights = embree::reduce_add(currentWeightsVec); - float previousWeights = embree::reduce_add(previousWeightsVec); - float sumPreviousWeights = embree::reduce_add(sumPreviousWeightsVec); - // float sumCurrentWeights = embree::reduce_add(sumCurrentWeightsVec); - // std::cout << "sumCurrentWeights = "<< sumCurrentWeights << "\t currentStats.sumWeights = "<< currentStats.sumWeights << std::endl; - // std::cout << "sumPreviousWeights = " << sumPreviousWeights << "\t previousStats.sumWeights = "<< previousStats.sumWeights << std::endl; - - embree::vfloat sumTmpVec = 0.f; - // calcualting the normalization factor for the weights of the updated components - float inv_currentWeights = (sumPreviousWeights - previousWeights) / currentWeights; - // std::cout << "sumPreviousWeights = " << sumPreviousWeights << "\t previousWeights = "<< previousWeights << "\t currentWeights = "<< currentWeights << "\t inv_currentWeights - // = "<< inv_currentWeights << std::endl; - for (int k = 0; k < cnt; k++) - { - currentStats.sumOfWeightedStats[k] = select(mask.mask[k], currentStats.sumOfWeightedStats[k] * inv_currentWeights, currentStats.sumOfWeightedStats[k]); - currentStats.sumOfWeightedDirections[k].x = select(mask.mask[k], currentStats.sumOfWeightedDirections[k].x * inv_currentWeights, currentStats.sumOfWeightedDirections[k].x); - currentStats.sumOfWeightedDirections[k].y = select(mask.mask[k], currentStats.sumOfWeightedDirections[k].y * inv_currentWeights, currentStats.sumOfWeightedDirections[k].y); - currentStats.sumOfWeightedDirections[k].z = select(mask.mask[k], currentStats.sumOfWeightedDirections[k].z * inv_currentWeights, currentStats.sumOfWeightedDirections[k].z); - sumTmpVec += currentStats.sumOfWeightedStats[k]; - - OPENPGL_ASSERT(embree::isvalid(currentStats.sumOfWeightedStats[k])); - OPENPGL_ASSERT(embree::isvalid(currentStats.sumOfWeightedDirections[k].x)); - OPENPGL_ASSERT(embree::isvalid(currentStats.sumOfWeightedDirections[k].y)); - OPENPGL_ASSERT(embree::isvalid(currentStats.sumOfWeightedDirections[k].z)); - } - - // since only some partial components are updated - // the overall stats are the same as the previous stats - currentStats.sumWeights = previousStats.sumWeights; - currentStats.numSamples = previousStats.numSamples; - currentStats.overallNumSamples = previousStats.overallNumSamples; - currentStats.norm = previousStats.norm; - currentStats.inv_norm = previousStats.inv_norm; -} -#endif template float ParallaxAwareVonMisesFisherWeightedEMFactory::weightedExpectationStep(VMM &vmm, SufficientStatistics &stats, UnassignedSamplesStatistics &unassignedStats,