|| || || || || ||

---

Chi2DistanceSSE.h

Go to the documentation of this file.

#ifndef Impala_Core_Vector_Chi2DistanceSSE_h
#define Impala_Core_Vector_Chi2DistanceSSE_h

#include "Core/Vector/VectorTem.h"
// for includes, defines
#include "Core/Matrix/MatNorm2DistSSE.h"

namespace Impala
{
namespace Core
{
namespace Vector
{


double
Chi2DistanceSSE(const VectorTem<double>& v1, const VectorTem<double>& v2)
{
    // please do not add ILOG_VAR to this function
    int length=v1.Size();
    double* v1Copy = 0;
    double* v2Copy = 0;
    const double* v1Data = v1.GetData();
    const double* v2Data = v2.GetData();
    if(!POINTER_ALIGNED(v1Data))
    {
        // make a copy
        v1Copy = new double[length];
        for(int i = 0; i < length; i++)
            v1Copy[i] = v1Data[i];
        v1Data = v1Copy;
    }
    if(!POINTER_ALIGNED(v2Data))
    {
        // make a copy
        v2Copy = new double[length];
        for(int i = 0; i < length; i++)
            v2Copy[i] = v2Data[i];
        v2Data = v2Copy;
    }
    __m128d* baseA = (__m128d*)v1Data;
    __m128d* baseB = (__m128d*)v2Data;
    const int SSELength = IntAlignDown(length, 2) / 2;
    const __m128d zero = _mm_setzero_pd();
    __m128d accu = _mm_setzero_pd();
    for (unsigned int k = 0; k < SSELength; k++)
    {
        __m128d As = baseA[k];
        __m128d Bs = baseB[k];
        __m128d tmp = _mm_sub_pd(As, Bs);
        // sum += (As - Bs) * (As - Bs) / (As + Bs);
        tmp =_mm_div_pd(_mm_mul_pd(tmp, tmp), _mm_add_pd(As, Bs));
        // fix the NaN for 0/0
        tmp = _mm_max_pd(tmp, zero);
        accu = _mm_add_pd(accu, tmp);
    }
    // add the components of sum together
    __m128d shuffle = _mm_shuffle_pd(accu, accu, _MM_SHUFFLE2(0, 1));
    __m128d result  = _mm_add_pd(accu, shuffle) ;
    double accumulator = _mm_cvtsd_f64(result);

    for(int i=SSELength*2 ; i<length ; i++)
    {
        double sum = v1[i] + v2[i];
        if(sum)
            accumulator += (v1[i] - v2[i]) * (v1[i] - v2[i]) / sum;
    }
    if(v1Copy) delete[] v1Copy;
    if(v2Copy) delete[] v2Copy;

    // make it between 0..1, instead of 0..2 (this statement is only true if the input vectors sum to one)
    return accumulator / 2.0;
}


} // namespace Vector
} // namespace Core
} // namespace Impala


#endif

---