vertical_segment.cpp 4.79 KB
Newer Older
1 2
/* ALTA --- Analysis of Bidirectional Reflectance Distribution Functions

3
   Copyright (C) 2013, 2014, 2015, 2016, 2017 Inria
4 5 6 7 8 9 10

   This file is part of ALTA.

   This Source Code Form is subject to the terms of the Mozilla Public
   License, v. 2.0.  If a copy of the MPL was not distributed with this
   file, You can obtain one at http://mozilla.org/MPL/2.0/.  */

11
#include "vertical_segment.h"
12
#include "data_storage.h"
13 14 15 16 17 18 19

#include <string>
#include <sstream>
#include <iostream>
#include <fstream>
#include <algorithm>
#include <cmath>
20
#include <cassert>
21
#include <cstring>
22

23
using namespace alta;
24
using namespace Eigen;
25

26
//#define RELATIVE_ERROR
27

28 29
// Note: For some reason returning an rvalue here doesn't work, so let's hope
// RVO comes into play.
30
static vec data_min(Ref<MatrixXd> data)
31
{
32 33
    vec min = VectorXd::Constant(data.rows(),
                                 std::numeric_limits<double>::max());
34

35
    for(int i = 0; i < data.cols(); i++)
36
    {
37 38
        auto col = data.col(i);
        for (int k = 0; k < data.rows(); ++k)
39
        {
40
            min[k] = std::min(min[k], col[k]);
41 42 43 44
        }
    }

    return min;
45 46
}

47
static vec data_max(Ref<MatrixXd> data)
48
{
49 50
    vec max = VectorXd::Constant(data.rows(),
                                 -std::numeric_limits<double>::max());
51

52
    for(int i = 0; i < data.cols(); i++)
53
    {
54 55
        auto col = data.col(i);
        for (int k = 0; k < data.rows(); ++k)
56
        {
57
            max[k] = std::max(max[k], col[k]);
58 59 60 61 62 63
        }
    }

    return max;
}

64 65 66
// Return a matrix view of PTR showing only the 'dimX' first rows of that
// matrix.
static Ref<MatrixXd>
67 68
x_view(double *ptr, size_t cols, const parameters& params,
       vertical_segment::ci_kind kind)
69
{
70 71 72
    auto stride = params.dimX() + params.dimY()
        + vertical_segment::confidence_interval_columns(kind, params);

73
    return Map<MatrixXd, 0, OuterStride<> >(ptr, params.dimX(), cols,
74
                                            OuterStride<>(stride));
75
}
76

77
vertical_segment::vertical_segment(const parameters& params,
78
                                   size_t size,
79 80
                                   std::shared_ptr<double> input_data,
                                   ci_kind kind)
81
    : data(params, size,
82 83 84 85 86
           data_min(x_view(input_data.get(), size, params, kind)),
           data_max(x_view(input_data.get(), size, params, kind))),
      _data(input_data),
      _ci_kind(kind),
      _is_absolute(true), _dt(0.1)
87 88 89
{
}

90 91 92 93 94 95
// Work around the lack of array support in C++11's 'shared_ptr'.
static void delete_array(double *thing)
{
    delete[] thing;
}

96
vertical_segment::vertical_segment(const parameters& params, unsigned int rows):
97 98 99 100 101
    // Create a data object with a "value-initialized" (i.e., zeroed) array.
    vertical_segment(params, rows,
                     std::shared_ptr<double>
                     (new double[rows * (params.dimX() + 3 * params.dimY())]{},
                      delete_array))
102 103 104
{
}

105

106
void vertical_segment::get(int i, vec& x, vec& yl, vec& yu) const
107
{
108 109 110
    // Make sure we have the lower and upper bounds of Y.
    assert(confidence_interval_kind() == ASYMMETRICAL_CONFIDENCE_INTERVAL);

111
    auto row = matrix_view().row(i);
112

113
    x = row.segment(0, _parameters.dimX());
114

115
    auto y = row.segment(_parameters.dimX(), _parameters.dimY());
116 117

    yl = row.segment(_parameters.dimX() + _parameters.dimY(),
118
                     _parameters.dimY());
119 120

    yu = row.segment(_parameters.dimX() + 2 * _parameters.dimY(),
121
                     _parameters.dimY());
122
}
123

124 125
void vertical_segment::get(int i, vec& yl, vec& yu) const
{
126 127
    vec x;
    get(i, x, yl, yu);
128
}
129

130
vec vertical_segment::get(int i) const
131
{
132
    return data_view().row(i);
133 134
}

135 136
void vertical_segment::set(int i, const vec& x)
{
137 138
   // Check if the input data 'x' has the size of a vertical segment (i.e. dimX+3*dimY),
   // if it only has the size of a single value, then create the segment.
139 140 141
   if(x.size() == column_number()
      || x.size() == _parameters.dimX() + _parameters.dimY())
   {
142
       auto row = matrix_view().row(i);
143
       row.head(x.size()) = x;
144 145 146 147
   } else {
      std::cerr << "<<ERROR>> Passing an incorrect element to vertical_segment::set" << std::endl;
      throw;
   }
148 149
}

150
vec vertical_segment::vs(const vec& x) const {
151
   vec y(_parameters.dimX() + 3*_parameters.dimY());
152 153

   // Copy the head of each vector
154
   y.head(_parameters.dimX() + _parameters.dimY()) = x.head(_parameters.dimX() + _parameters.dimY());
155
   
156 157
   for(unsigned int i=0; i<_parameters.dimY(); ++i) {
      const double val = x[i + _parameters.dimX()];
158
      if(_is_absolute) {
159 160
         y[i + _parameters.dimX()+1*_parameters.dimY()] = val - _dt;
         y[i + _parameters.dimX()+2*_parameters.dimY()] = val + _dt;
161
      } else {
162 163
         y[i + _parameters.dimX()+1*_parameters.dimY()] = val * (1.0 - _dt);
         y[i + _parameters.dimX()+2*_parameters.dimY()] = val * (1.0 + _dt);
164
      }
165 166 167 168
   }

   return y;
}