// This file is part of HDRip.
//
// HDRip is free software: you can redistribute it and/or modify it
// under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// HDRip is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with HDRip. If not, see .
//
// HDRip project
// Author : Rémi Synave
// Contact : remi.synave@univ-littoral.fr
#include "pch.h"
#include
#include
#include "Utils.hpp"
#include "MT_interpolation.hpp"
#ifdef _MT_
void* interp_MT(void* arg)
{
MT_interpolation* a = (MT_interpolation*)arg;
const float* xData = a->xData;
float* xInterp = a->xInterp;
for (unsigned int i = 0; i < a->length; i++)
{
unsigned int index = Utils::search(xData[i], a->xd, 0, a->length_xd_yd - 1);
xInterp[i] = Utils::interp(xData[i], a->xd[index], a->yd[index], a->xd[index + 1], a->yd[index + 1]);
}
return arg;
}
float* Utils::interp(const float* x, const unsigned int tx, const Eigen::RowVectorXf& xd, const Eigen::RowVectorXf& yd)
{
float* xInterp = new float[tx];
std::thread tab_t[_MT_];
MT_interpolation tab_a[_MT_];
unsigned int id;
unsigned int tab_length = tx;
unsigned int block_size = tab_length / _MT_;
for (id = 0; id < _MT_; id++) {
tab_a[id].xInterp = xInterp + (id * block_size);
tab_a[id].xData = x + (id * block_size);
tab_a[id].length = block_size;
tab_a[id].xd = xd.data();
tab_a[id].yd = yd.data();
tab_a[id].length_xd_yd = (unsigned int)(xd.size());
if (id == (_MT_ - 1))
tab_a[id].length = tab_length - ((_MT_ - 1) * block_size);
tab_t[id] = std::thread(interp_MT, (void*)(tab_a + id));
}
for (id = 0; id < _MT_; id++) {
tab_t[id].join();
}
return xInterp;
}
#else
float* Utils::interp(const float* x, const unsigned int tx, const Eigen::RowVectorXf& xd, const Eigen::RowVectorXf& yd)
{
float* xInterp = new float[tx];
for (unsigned int i = 0; i < tx; i++)
{
unsigned int index = Utils::search(x[i], xd);
xInterp[i] = interp(x[i], xd(index), yd(index), xd(index + 1), yd(index + 1));
}
return xInterp;
}
#endif
float Utils::interp(float x, float xd1, float yd1, float xd2, float yd2)
{
float delta_xd2_xd1 = xd2 - xd1;
float delta_x_xd1 = x - xd1;
float delta_yd2_yd1 = yd2 - yd1;
return yd1 + delta_yd2_yd1 * (delta_x_xd1 / delta_xd2_xd1);
}
// Maybe this method is useless and the index can be determined in the interp metho.
unsigned int Utils::search(float to_search, const Eigen::RowVectorXf& data)
{
// In case of approximation
if (to_search <= data[0])
return 0;
if (to_search >= data[data.size() - 1])
return (unsigned int)(data.size() - 2);
unsigned int first = 0;
unsigned int last = (unsigned int)(data.size() - 1);
unsigned int index = (first + last) / 2;
while (!(to_search >= data[index] && to_search < data[index + 1]))
{
if (to_search < data[index])
last = index;
else
first = index;
index = (first + last) / 2;
}
return index;
}
// Maybe this method is useless and the index can be determined in the interp metho.
unsigned int Utils::search(float to_search, const float* data, unsigned int first, unsigned int last)
{
// In case of approximation
if (to_search <= data[0])
return 0;
if (to_search >= data[last])
return last - 1;
unsigned int index = (first + last) / 2;
while (!(to_search >= data[index] && to_search < data[index + 1]))
{
if (to_search < data[index])
last = index;
else
first = index;
index = (first + last) / 2;
}
return index;
}
// Parallelization is unless for this method because it is used in an already parallelized method.
float* Utils::NPlinearWeightMask(float* x, unsigned int length, float xMin, float xMax, float xTolerance)
{
float* y = new float[length];
for(unsigned int i=0;i(xMin-xTolerance)) && (x[i]<=xMin))
y[i]=(x[i]-(xMin-xTolerance)) / xTolerance;
for(unsigned int i=0;ixMin) && (x[i]<=xMax))
y[i]=1;
for(unsigned int i=0;ixMax) && (x[i]<=(xMax+xTolerance)))
y[i]=1-((x[i]-xMax) / xTolerance);
for(unsigned int i=0;i(xMax+xTolerance))
y[i]=0;
return y;
}
//#endif