← Back to branch summary

~valavanisalex/ubuntu/maverick/scidavis/fix-604811

« back to all changes in this revision

Viewing changes to scidavis/src/Interpolation.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Ruben Molina
  • Date: 2009-09-06 11:34:04 UTC
  • Revision ID: james.westby@ubuntu.com-20090906113404-4awaey82l3686w4q
Tags: upstream-0.2.3
Import upstream version 0.2.3

Show diffs side-by-side

added added

removed removed

Lines of Context:
scidavis/src/Interpolation.cpp
 
1
/***************************************************************************
 
2
    File                 : Interpolation.cpp
 
3
    Project              : SciDAVis
 
4
    --------------------------------------------------------------------
 
5
    Copyright            : (C) 2007 by Ion Vasilief
 
6
    Email (use @ for *)  : ion_vasilief*yahoo.fr
 
7
    Description          : Numerical interpolation of data sets
 
8
 
 
9
 ***************************************************************************/
 
10
 
 
11
/***************************************************************************
 
12
 *                                                                         *
 
13
 *  This program is free software; you can redistribute it and/or modify   *
 
14
 *  it under the terms of the GNU General Public License as published by   *
 
15
 *  the Free Software Foundation; either version 2 of the License, or      *
 
16
 *  (at your option) any later version.                                    *
 
17
 *                                                                         *
 
18
 *  This program is distributed in the hope that it will be useful,        *
 
19
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of         *
 
20
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the          *
 
21
 *  GNU General Public License for more details.                           *
 
22
 *                                                                         *
 
23
 *   You should have received a copy of the GNU General Public License     *
 
24
 *   along with this program; if not, write to the Free Software           *
 
25
 *   Foundation, Inc., 51 Franklin Street, Fifth Floor,                    *
 
26
 *   Boston, MA  02110-1301  USA                                           *
 
27
 *                                                                         *
 
28
 ***************************************************************************/
 
29
#include "Interpolation.h"
 
30
 
 
31
#include <QMessageBox>
 
32
 
 
33
#include <gsl/gsl_sort.h>
 
34
#include <gsl/gsl_spline.h>
 
35
#include <gsl/gsl_interp.h>
 
36
 
 
37
Interpolation::Interpolation(ApplicationWindow *parent, Graph *g, const QString& curveTitle, int m)
 
38
: Filter(parent, g)
 
39
{
 
40
        init(m);
 
41
        setDataFromCurve(curveTitle);
 
42
}
 
43
 
 
44
Interpolation::Interpolation(ApplicationWindow *parent, Graph *g, const QString& curveTitle,
 
45
                             double start, double end, int m)
 
46
: Filter(parent, g)
 
47
{
 
48
        init(m);
 
49
        setDataFromCurve(curveTitle, start, end);
 
50
}
 
51
 
 
52
void Interpolation::init(int m)
 
53
{
 
54
    if (m < 0 || m > 2)
 
55
    {
 
56
        QMessageBox::critical((ApplicationWindow *)parent(), tr("SciDAVis") + " - " + tr("Error"),
 
57
        tr("Unknown interpolation method. Valid values are: 0 - Linear, 1 - Cubic, 2 - Akima."));
 
58
        d_init_err = true;
 
59
        return;
 
60
    }
 
61
    d_method = m;
 
62
        switch(d_method)
 
63
        {
 
64
                case 0:
 
65
                        setName(tr("Linear") + tr("Int"));
 
66
                        d_explanation = tr("Linear") + " " + tr("Interpolation");
 
67
                        break;
 
68
                case 1:
 
69
                        setName(tr("Cubic") + tr("Int"));
 
70
                        d_explanation = tr("Cubic") + " " + tr("Interpolation");
 
71
                        break;
 
72
                case 2:
 
73
                        setName(tr("Akima") + tr("Int"));
 
74
                        d_explanation = tr("Akima") + " " + tr("Interpolation");
 
75
                        break;
 
76
        }
 
77
    d_sort_data = true;
 
78
    d_min_points = d_method + 3;
 
79
}
 
80
 
 
81
 
 
82
void Interpolation::setMethod(int m)
 
83
{
 
84
if (m < 0 || m > 2)
 
85
    {
 
86
        QMessageBox::critical((ApplicationWindow *)parent(), tr("Error"),
 
87
        tr("Unknown interpolation method, valid values are: 0 - Linear, 1 - Cubic, 2 - Akima."));
 
88
        d_init_err = true;
 
89
        return;
 
90
    }
 
91
int min_points = m + 3;
 
92
if (d_n < min_points)
 
93
        {
 
94
                QMessageBox::critical((ApplicationWindow *)parent(), tr("SciDAVis") + " - " + tr("Error"),
 
95
                                tr("You need at least %1 points in order to perform this operation!").arg(min_points));
 
96
        d_init_err = true;
 
97
        return;
 
98
        }
 
99
d_method = m;
 
100
d_min_points = min_points;
 
101
}
 
102
 
 
103
void Interpolation::calculateOutputData(double *x, double *y)
 
104
{
 
105
        gsl_interp_accel *acc = gsl_interp_accel_alloc ();
 
106
        const gsl_interp_type *method;
 
107
        switch(d_method)
 
108
        {
 
109
                case 0:
 
110
                        method = gsl_interp_linear;
 
111
                        break;
 
112
                case 1:
 
113
                        method = gsl_interp_cspline;
 
114
                        break;
 
115
                case 2:
 
116
                        method = gsl_interp_akima;
 
117
                        break;
 
118
        }
 
119
 
 
120
        gsl_spline *interp = gsl_spline_alloc (method, d_n);
 
121
        gsl_spline_init (interp, d_x, d_y, d_n);
 
122
 
 
123
    double step = (d_to - d_from)/(double)(d_points - 1);
 
124
    for (int j = 0; j < d_points; j++)
 
125
        {
 
126
           x[j] = d_from + j*step;
 
127
           y[j] = gsl_spline_eval (interp, x[j], acc);
 
128
        }
 
129
 
 
130
        gsl_spline_free (interp);
 
131
        gsl_interp_accel_free (acc);
 
132
}
 
133
 
 
134
int Interpolation::sortedCurveData(QwtPlotCurve *c, double start, double end, double **x, double **y)
 
135
{
 
136
    if (!c || c->rtti() != QwtPlotItem::Rtti_PlotCurve)
 
137
        return 0;
 
138
 
 
139
    int i_start = 0, i_end = c->dataSize();
 
140
    for (int i = 0; i < i_end; i++)
 
141
            if (c->x(i) > start && i)
 
142
        {
 
143
              i_start = i - 1;
 
144
          break;
 
145
        }
 
146
    for (int i = i_end-1; i >= 0; i--)
 
147
            if (c->x(i) < end && i < c->dataSize())
 
148
        {
 
149
              i_end = i + 1;
 
150
          break;
 
151
        }
 
152
    int n = i_end - i_start + 1;
 
153
    (*x) = new double[n];
 
154
    (*y) = new double[n];
 
155
    double *xtemp = new double[n];
 
156
    double *ytemp = new double[n];
 
157
 
 
158
        double pr_x;
 
159
        int j=0;
 
160
    for (int i = i_start; i <= i_end; i++)
 
161
    {
 
162
        xtemp[j] = c->x(i);
 
163
        if (xtemp[j] == pr_x)
 
164
        {
 
165
            delete (*x);
 
166
            delete (*y);
 
167
            return -1;//this kind of data causes division by zero in GSL interpolation routines
 
168
        }
 
169
        pr_x = xtemp[j];
 
170
        ytemp[j++] = c->y(i);
 
171
    }
 
172
    size_t *p = new size_t[n];
 
173
    gsl_sort_index(p, xtemp, 1, n);
 
174
    for (int i=0; i<n; i++)
 
175
    {
 
176
        (*x)[i] = xtemp[p[i]];
 
177
            (*y)[i] = ytemp[p[i]];
 
178
    }
 
179
    delete[] xtemp;
 
180
    delete[] ytemp;
 
181
    delete[] p;
 
182
    return n;
 
183
}