GaussianOffsetsCalculator

GaussianOffsetsCalculator.cs

using System;
using System.Collections.Generic;
using System.Linq;

namespace GaussianOffsetsCalculator
{
    class Program
    {

        private static long[] _fracCache = new long[30];

        private static void InitFracCache()
        {
            long frac = 1;
            _fracCache[0] = 1;
            for (int i = 1; i < _fracCache.Length; ++i)
            {
                frac *= i;
                _fracCache[i] = frac;
            }
        }

        //TODO: 看起来这里不能用非递推算法，太容易溢出了……
        private static long C(int n, int k)
        {
            return _fracCache[n] / (_fracCache[k] * _fracCache[n - k]);
        }

        private static long[] GetRowNOfPascalTriangle(int n)
        {
            var result = new long[n + 1];
            for (int i = 0; i <= n; ++i)
            {
                result[i] = C(n, i);
            }
            return result;
        }

        private static float[] RemoveMargin(float[] weights, int margin)
        {
            return weights.Skip(margin).Take(weights.Length - 2 * margin).ToArray();
        }
        private static float[] TryGetRawGaussianKernelWeights(int n)
        {
            float[] result = new float[n + 1];
            long[] row = GetRowNOfPascalTriangle(n);
            int coefficient = 1 << n;
            for (int i = 0; i <= n; ++i)
            {
                result[i] = (float)row[i] / (float)coefficient;
            }
            return result;
        }

        private static float[] GetRawGaussianKernelWeights(int n, int expected)
        {
            var weights = TryGetRawGaussianKernelWeights(n);
            return RemoveMargin(weights, (n + 1 - expected) / 2);
        }

        private static (float[] optedOffsets, float[] optedWeights) OptimizeWeights(float[] weights)
        {
            int start = (weights.Length + 1) / 2;
            int outputLength = 1 + ((weights.Length - 1) / 2 + 1) / 2;
            int length = (weights.Length - start ) / 2;
            float[] optedOffsets = new float[outputLength];
            float[] optedWeights = new float[outputLength];
            
            optedOffsets[0] = 0.0f;
            optedWeights[0] = weights[start - 1];
            for (int i = 0; i < length * 2; i += 2)
            {
                float weight1 = weights[start + i];
                float weight2 = weights[start + i + 1];
                float offset1 = i + 1;
                float offset2 = i + 2;

                optedWeights[i / 2 + 1] = weight1 + weight2;
                optedOffsets[i / 2 + 1] = (offset1 * weight1 + offset2 * weight2) / (weight1 + weight2);
            }

            if (outputLength > length + 1)
            {
                optedWeights[outputLength - 1] = weights[weights.Length - 1];
                optedOffsets[outputLength - 1] = (weights.Length - 1) / 2;
            }
            return (optedOffsets, optedWeights);
        }

        static void Main(string[] args)
        {
            InitFracCache();
            // long[] row = GetRowNOfPascalTriangle(14);
            // Console.WriteLine(string.Join(", ", row));

            float[] rawWeights = GetRawGaussianKernelWeights(12, 11);
            Console.WriteLine(string.Join(", ", rawWeights));

            var (optedOffsets, optedWeights) = OptimizeWeights(rawWeights);
            Console.WriteLine($"weights: {string.Join(", ", optedWeights)}");
            Console.WriteLine($"offsets: {string.Join(", ", optedOffsets)}");
        }
    }
}