diff --git a/src/backend/core/AbstractColumn.cpp b/src/backend/core/AbstractColumn.cpp index 0aebf8c85..e96e23e04 100644 --- a/src/backend/core/AbstractColumn.cpp +++ b/src/backend/core/AbstractColumn.cpp @@ -1,708 +1,710 @@ /*************************************************************************** File : AbstractColumn.cpp Project : LabPlot Description : Interface definition for data with column logic -------------------------------------------------------------------- Copyright : (C) 2007,2008 Tilman Benkert (thzs@gmx.net) Copyright : (C) 2017 Stefan Gerlach (stefan.gerlach@uni.kn) ***************************************************************************/ /*************************************************************************** * * * This program 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 2 of the License, or * * (at your option) any later version. * * * * This program 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 this program; if not, write to the Free Software * * Foundation, Inc., 51 Franklin Street, Fifth Floor, * * Boston, MA 02110-1301 USA * * * ***************************************************************************/ #include "backend/core/AbstractColumn.h" #include "backend/core/AbstractColumnPrivate.h" #include "backend/core/abstractcolumncommands.h" #include "backend/lib/XmlStreamReader.h" #include "backend/lib/SignallingUndoCommand.h" #include #include #include #include /** * \class AbstractColumn * \brief Interface definition for data with column logic * * This is an abstract base class for column-based data, * i.e. mathematically a vector or technically a 1D array or list. * It only defines the interface but has no data members itself. * * Classes derived from this are typically table columns or outputs * of filters which can be chained between table columns and plots. * From the point of view of the plot functions there will be no difference * between a table column and a filter output since both use this interface. * * Classes derived from this will either store a * vector with entries of one certain data type, e.g. double, QString, * QDateTime, or generate such values on demand. To determine the data * type of a class derived from this, use the columnMode() function. * AbstractColumn defines all access functions for all supported data * types but only those corresponding to the return value of columnMode() * will return a meaningful value. Calling functions not belonging to * the data type of the column is safe, but will do nothing (writing * function) or return some default value (reading functions). * * This class also defines all signals which indicate a data change. * Any class whose output values are subject to change over time must emit * the according signals. These signals notify any object working with the * column before and after a change of the column. * In some cases it will be necessary for a class using * the column to connect aboutToBeDestroyed(), to react * to a column's deletion, e.g. a filter's reaction to a * table deletion. * * All writing functions have a "do nothing" standard implementation to * make deriving a read-only class very easy without bothering about the * writing interface. */ /** * \brief Ctor * * \param name the column name (= aspect name) */ AbstractColumn::AbstractColumn(const QString &name) : AbstractAspect(name), m_abstract_column_private( new AbstractColumnPrivate(this) ) { } AbstractColumn::~AbstractColumn() { emit aboutToBeDestroyed(this); delete m_abstract_column_private; } QStringList AbstractColumn::dateFormats() { QStringList dates; dates << "yyyy-MM-dd"; dates << "yyyy/MM/dd"; dates << "dd/MM/yyyy"; dates << "dd/MM/yy"; dates << "dd.MM.yyyy"; dates << "dd.MM.yy"; dates << "MM/yyyy"; dates << "dd.MM."; dates << "yyyyMMdd"; return dates; } QStringList AbstractColumn::timeFormats() { QStringList times; times << "hh"; times << "hh ap"; times << "hh:mm"; times << "hh:mm ap"; times << "hh:mm:ss"; times << "hh:mm:ss.zzz"; times << "hh:mm:ss:zzz"; times << "mm:ss.zzz"; times << "hhmmss"; return times; } QStringList AbstractColumn::dateTimeFormats() { QStringList dateTimes; // any combination of date and times for (auto d: dateFormats()) dateTimes << d; for (auto t: timeFormats()) dateTimes << t; for (auto d: dateFormats()) for (auto t: timeFormats()) dateTimes << d + ' ' + t; return dateTimes; } /** * \brief Convenience method for mode-dependent icon */ QIcon AbstractColumn::iconForMode(ColumnMode mode) { switch (mode) { case AbstractColumn::Numeric: case AbstractColumn::Integer: break; case AbstractColumn::Text: return QIcon::fromTheme("draw-text"); case AbstractColumn::DateTime: case AbstractColumn::Month: case AbstractColumn::Day: return QIcon::fromTheme("chronometer"); } return QIcon::fromTheme("x-shape-text"); } /** * \fn bool AbstractColumn::isReadOnly() const * \brief Return whether the object is read-only */ /** * \fn AbstractColumn::ColumnMode AbstractColumn::columnMode() const * \brief Return the column mode * * This function is most used by tables but can also be used * by plots. The column mode specifies how to interpret * the values in the column additional to the data type. */ /** * \brief Set the column mode * * This sets the column mode and, if * necessary, converts it to another datatype. */ void AbstractColumn::setColumnMode(AbstractColumn::ColumnMode) {} /** * \brief Copy another column of the same type * * This function will return false if the data type * of 'other' is not the same as the type of 'this'. * Use a filter to convert a column to another type. */ bool AbstractColumn::copy(const AbstractColumn *other) { Q_UNUSED(other) return false; } /** * \brief Copies part of another column of the same type * * This function will return false if the data type * of 'other' is not the same as the type of 'this'. * \param source pointer to the column to copy * \param source_start first row to copy in the column to copy * \param destination_start first row to copy in * \param num_rows the number of rows to copy */ bool AbstractColumn::copy(const AbstractColumn *source, int source_start, int destination_start, int num_rows) { Q_UNUSED(source) Q_UNUSED(source_start) Q_UNUSED(destination_start) Q_UNUSED(num_rows) return false; } /** * \fn int AbstractColumn::rowCount() const * \brief Return the data vector size */ /** * \brief Insert some empty (or initialized with invalid values) rows */ void AbstractColumn::insertRows(int before, int count) { beginMacro( i18np("%1: insert 1 row", "%1: insert %2 rows", name(), count) ); exec(new SignallingUndoCommand("pre-signal", this, "rowsAboutToBeInserted", "rowsRemoved", Q_ARG(const AbstractColumn*,this), Q_ARG(int,before), Q_ARG(int,count))); handleRowInsertion(before, count); exec(new SignallingUndoCommand("post-signal", this, "rowsInserted", "rowsAboutToBeRemoved", Q_ARG(const AbstractColumn*,this), Q_ARG(int,before), Q_ARG(int,count))); endMacro(); } void AbstractColumn::handleRowInsertion(int before, int count) { exec(new AbstractColumnInsertRowsCmd(this, before, count)); } /** * \brief Remove 'count' rows starting from row 'first' */ void AbstractColumn::removeRows(int first, int count) { beginMacro( i18np("%1: remove 1 row", "%1: remove %2 rows", name(), count) ); exec(new SignallingUndoCommand("change signal", this, "rowsAboutToBeRemoved", "rowsInserted", Q_ARG(const AbstractColumn*,this), Q_ARG(int,first), Q_ARG(int,count))); handleRowRemoval(first, count); exec(new SignallingUndoCommand("change signal", this, "rowsRemoved", "rowsAboutToBeInserted", Q_ARG(const AbstractColumn*,this), Q_ARG(int,first), Q_ARG(int,count))); endMacro(); } void AbstractColumn::handleRowRemoval(int first, int count) { exec(new AbstractColumnRemoveRowsCmd(this, first, count)); } /** * \fn AbstractColumn::PlotDesignation AbstractColumn::plotDesignation() const * \brief Return the column plot designation */ /** * \brief Set the column plot designation */ void AbstractColumn::setPlotDesignation(AbstractColumn::PlotDesignation pd) { Q_UNUSED(pd) } /** * \brief Clear the whole column */ void AbstractColumn::clear() {} /** * \brief Convenience method for mode-independent testing of validity */ bool AbstractColumn::isValid(int row) const { switch (columnMode()) { case AbstractColumn::Numeric: return !std::isnan(valueAt(row)); case AbstractColumn::Integer: // there is no invalid integer return true; case AbstractColumn::Text: return !textAt(row).isNull(); case AbstractColumn::DateTime: case AbstractColumn::Month: case AbstractColumn::Day: return dateTimeAt(row).isValid(); } return false; } //////////////////////////////////////////////////////////////////////////////////////////////////// //! \name IntervalAttribute related functions //@{ //////////////////////////////////////////////////////////////////////////////////////////////////// /** * \brief Return whether a certain row is masked */ bool AbstractColumn::isMasked(int row) const { return m_abstract_column_private->m_masking.isSet(row); } /** * \brief Return whether a certain interval of rows is fully masked */ bool AbstractColumn::isMasked(Interval i) const { return m_abstract_column_private->m_masking.isSet(i); } /** * \brief Return all intervals of masked rows */ QList< Interval > AbstractColumn::maskedIntervals() const { return m_abstract_column_private->m_masking.intervals(); } /** * \brief Clear all masking information */ void AbstractColumn::clearMasks() { exec(new AbstractColumnClearMasksCmd(m_abstract_column_private), "maskingAboutToChange", "maskingChanged", Q_ARG(const AbstractColumn*,this)); } /** * \brief Set an interval masked * * \param i the interval * \param mask true: mask, false: unmask */ void AbstractColumn::setMasked(Interval i, bool mask) { exec(new AbstractColumnSetMaskedCmd(m_abstract_column_private, i, mask), "maskingAboutToChange", "maskingChanged", Q_ARG(const AbstractColumn*,this)); } /** * \brief Overloaded function for convenience */ void AbstractColumn::setMasked(int row, bool mask) { setMasked(Interval(row,row), mask); } //////////////////////////////////////////////////////////////////////////////////////////////////// //@} //////////////////////////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////////////////////////// //! \name Formula related functions //@{ //////////////////////////////////////////////////////////////////////////////////////////////////// /** * \brief Return the formula associated with row 'row' */ QString AbstractColumn::formula(int row) const { Q_UNUSED(row); return QString(); } /** * \brief Return the intervals that have associated formulas * * This can be used to make a list of formulas with their intervals. * Here is some example code: * * \code * QStringList list; * QList< Interval > intervals = my_column.formulaIntervals(); * foreach(Interval interval, intervals) * list << QString(interval.toString() + ": " + my_column.formula(interval.start())); * \endcode */ QList< Interval > AbstractColumn::formulaIntervals() const { return QList< Interval >(); } /** * \brief Set a formula string for an interval of rows */ void AbstractColumn::setFormula(Interval i, QString formula) { Q_UNUSED(i) Q_UNUSED(formula) } /** * \brief Overloaded function for convenience */ void AbstractColumn::setFormula(int row, QString formula) { Q_UNUSED(row) Q_UNUSED(formula) } /** * \brief Clear all formulas */ void AbstractColumn::clearFormulas() {}; //////////////////////////////////////////////////////////////////////////////////////////////////// //@} //////////////////////////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////////////////////////// //! \name type specific functions //@{ //////////////////////////////////////////////////////////////////////////////////////////////////// /** * \brief Return the content of row 'row'. * * Use this only when columnMode() is Text */ QString AbstractColumn::textAt(int row) const { Q_UNUSED(row); return ""; } /** * \brief Set the content of row 'row' * * Use this only when columnMode() is Text */ void AbstractColumn::setTextAt(int row, const QString& new_value) { Q_UNUSED(row) Q_UNUSED(new_value) } /** * \brief Replace a range of values * * Use this only when columnMode() is Text */ void AbstractColumn::replaceTexts(int first, const QVector& new_values) { Q_UNUSED(first) Q_UNUSED(new_values) }; /** * \brief Return the date part of row 'row' * * Use this only when columnMode() is DateTime, Month or Day */ QDate AbstractColumn::dateAt(int row) const { Q_UNUSED(row); return QDate(); } /** * \brief Set the content of row 'row' * * Use this only when columnMode() is DateTime, Month or Day */ void AbstractColumn::setDateAt(int row, const QDate& new_value) { Q_UNUSED(row) Q_UNUSED(new_value) }; /** * \brief Return the time part of row 'row' * * Use this only when columnMode() is DateTime, Month or Day */ QTime AbstractColumn::timeAt(int row) const { Q_UNUSED(row); return QTime(); } /** * \brief Set the content of row 'row' * * Use this only when columnMode() is DateTime, Month or Day */ void AbstractColumn::setTimeAt(int row, const QTime& new_value) { Q_UNUSED(row) Q_UNUSED(new_value) } /** * \brief Return the QDateTime in row 'row' * * Use this only when columnMode() is DateTime, Month or Day */ QDateTime AbstractColumn::dateTimeAt(int row) const { Q_UNUSED(row); return QDateTime(); } /** * \brief Set the content of row 'row' * * Use this only when columnMode() is DateTime, Month or Day */ void AbstractColumn::setDateTimeAt(int row, const QDateTime& new_value) { Q_UNUSED(row) Q_UNUSED(new_value) }; /** * \brief Replace a range of values * * Use this only when columnMode() is DateTime, Month or Day */ void AbstractColumn::replaceDateTimes(int first, const QVector& new_values) { Q_UNUSED(first) Q_UNUSED(new_values) }; /** * \brief Return the double value in row 'row' * * Use this only when columnMode() is Numeric */ double AbstractColumn::valueAt(int row) const { Q_UNUSED(row); return NAN; } /** * \brief Set the content of row 'row' * * Use this only when columnMode() is Numeric */ void AbstractColumn::setValueAt(int row, double new_value) { Q_UNUSED(row) Q_UNUSED(new_value) }; /** * \brief Replace a range of values * * Use this only when columnMode() is Numeric */ void AbstractColumn::replaceValues(int first, const QVector& new_values) { Q_UNUSED(first) Q_UNUSED(new_values) } /** * \brief Return the integer value in row 'row' * * Use this only when columnMode() is Integer */ int AbstractColumn::integerAt(int row) const { Q_UNUSED(row); return 42; } /** * \brief Set the content of row 'row' * * Use this only when columnMode() is Integer */ void AbstractColumn::setIntegerAt(int row, int new_value) { Q_UNUSED(row) Q_UNUSED(new_value) }; /** * \brief Replace a range of values * * Use this only when columnMode() is Integer */ void AbstractColumn::replaceInteger(int first, const QVector& new_values) { Q_UNUSED(first) Q_UNUSED(new_values) } /**********************************************************************/ double AbstractColumn::minimum(int count) const { + Q_UNUSED(count); return -INFINITY; } double AbstractColumn::maximum(int count) const { + Q_UNUSED(count); return INFINITY; } //////////////////////////////////////////////////////////////////////////////////////////////////// //@} //////////////////////////////////////////////////////////////////////////////////////////////////// /** * \fn void AbstractColumn::plotDesignationAboutToChange(const AbstractColumn *source) * \brief Column plot designation will be changed * * 'source' is always the this pointer of the column that * emitted this signal. This way it's easier to use * one handler for lots of columns. */ /** * \fn void AbstractColumn::plotDesignationChanged(const AbstractColumn *source) * \brief Column plot designation changed * * 'source' is always the this pointer of the column that * emitted this signal. This way it's easier to use * one handler for lots of columns. */ /** * \fn void AbstractColumn::modeAboutToChange(const AbstractColumn *source) * \brief Column mode (possibly also the data type) will be changed * * 'source' is always the this pointer of the column that * emitted this signal. This way it's easier to use * one handler for lots of columns. */ /** * \fn void AbstractColumn::modeChanged(const AbstractColumn *source) * \brief Column mode (possibly also the data type) changed * * 'source' is always the this pointer of the column that * emitted this signal. This way it's easier to use * one handler for lots of columns. */ /** * \fn void AbstractColumn::dataAboutToChange(const AbstractColumn *source) * \brief Data of the column will be changed * * 'source' is always the this pointer of the column that * emitted this signal. This way it's easier to use * one handler for lots of columns. */ /** * \fn void AbstractColumn::dataChanged(const AbstractColumn *source) * \brief Data of the column has changed * * Important: When data has changed also the number * of rows in the column may have changed without * any other signal emission. * 'source' is always the this pointer of the column that * emitted this signal. This way it's easier to use * one handler for lots of columns. */ /** * \fn void AbstractColumn::rowsAboutToBeInserted(const AbstractColumn *source, int before, int count) * \brief Rows will be inserted * * \param source the column that emitted the signal * \param before the row to insert before * \param count the number of rows to be inserted */ /** * \fn void AbstractColumn::rowsInserted(const AbstractColumn *source, int before, int count) * \brief Rows have been inserted * * \param source the column that emitted the signal * \param before the row to insert before * \param count the number of rows to be inserted */ /** * \fn void AbstractColumn::rowsAboutToBeRemoved(const AbstractColumn *source, int first, int count) * \brief Rows will be deleted * * \param source the column that emitted the signal * \param first the first row to be deleted * \param count the number of rows to be deleted */ /** * \fn void AbstractColumn::rowsRemoved(const AbstractColumn *source, int first, int count) * \brief Rows have been deleted * * \param source the column that emitted the signal * \param first the first row that was deleted * \param count the number of deleted rows */ /** * \fn void AbstractColumn::maskingAboutToChange(const AbstractColumn *source) * \brief Rows are about to be masked or unmasked */ /** * \fn void AbstractColumn::maskingChanged(const AbstractColumn *source) * \brief Rows have been masked or unmasked */ /** * \fn void AbstractColumn::aboutToBeDestroyed(const AbstractColumn *source) * \brief Emitted shortl before this column is deleted * * \param source the object emitting this signal * * This is needed by AbstractFilter. */ /** * \brief Read XML mask element */ bool AbstractColumn::XmlReadMask(XmlStreamReader *reader) { Q_ASSERT(reader->isStartElement() && reader->name() == "mask"); bool ok1, ok2; int start, end; start = reader->readAttributeInt("start_row", &ok1); end = reader->readAttributeInt("end_row", &ok2); if(!ok1 || !ok2) { reader->raiseError(i18n("invalid or missing start or end row")); return false; } setMasked(Interval(start,end)); if (!reader->skipToEndElement()) return false; return true; } /** * \brief Write XML mask element */ void AbstractColumn::XmlWriteMask(QXmlStreamWriter *writer) const { for (const auto& interval: maskedIntervals()) { writer->writeStartElement("mask"); writer->writeAttribute("start_row", QString::number(interval.start())); writer->writeAttribute("end_row", QString::number(interval.end())); writer->writeEndElement(); } } diff --git a/src/backend/nsl/nsl_geom_linesim.c b/src/backend/nsl/nsl_geom_linesim.c index 510b2b8ee..89b226448 100644 --- a/src/backend/nsl/nsl_geom_linesim.c +++ b/src/backend/nsl/nsl_geom_linesim.c @@ -1,618 +1,630 @@ /*************************************************************************** File : nsl_geom_linesim.c Project : LabPlot Description : NSL geometry line simplification functions -------------------------------------------------------------------- Copyright : (C) 2016 by Stefan Gerlach (stefan.gerlach@uni.kn) ***************************************************************************/ /*************************************************************************** * * * This program 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 2 of the License, or * * (at your option) any later version. * * * * This program 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 this program; if not, write to the Free Software * * Foundation, Inc., 51 Franklin Street, Fifth Floor, * * Boston, MA 02110-1301 USA * * * ***************************************************************************/ #include "nsl_geom_linesim.h" #include "nsl_geom.h" #include "nsl_common.h" #include "nsl_sort.h" #include "nsl_stats.h" const char* nsl_geom_linesim_type_name[] = {i18n("Douglas-Peucker (number)"), i18n("Douglas-Peucker (tolerance)"), i18n("Visvalingam-Whyatt"), i18n("Reumann-Witkam"), i18n("perpendicular distance"), i18n("n-th point"), i18n("radial distance"), i18n("Interpolation"), i18n("Opheim"), i18n("Lang")}; /*********** error calculation functions *********/ double nsl_geom_linesim_positional_error(const double xdata[], const double ydata[], const size_t n, const size_t index[]) { double dist = 0; size_t i = 0, j; /* i: index of index[] */ do { /*for every point not in index[] calculate distance to line*/ /*printf("i=%d (index[i]-index[i+1]=%d-%d)\n", i , index[i], index[i+1]);*/ for (j = 1; j < index[i+1]-index[i]; j++) { /*printf("i=%d: j=%d\n", i, j);*/ dist += nsl_geom_point_line_dist(xdata[index[i]], ydata[index[i]], xdata[index[i+1]], ydata[index[i+1]], xdata[index[i]+j], ydata[index[i]+j]); /*printf("dist = %g\n", dist);*/ } i++; } while(index[i] != n-1); return dist/(double)n; } double nsl_geom_linesim_positional_squared_error(const double xdata[], const double ydata[], const size_t n, const size_t index[]) { double dist = 0; size_t i = 0, j; /* i: index of index[] */ do { /*for every point not in index[] calculate distance to line*/ /*printf("i=%d (index[i]-index[i+1]=%d-%d)\n", i , index[i], index[i+1]);*/ for (j = 1; j < index[i+1]-index[i]; j++) { /*printf("i=%d: j=%d\n", i, j);*/ dist += pow(nsl_geom_point_line_dist(xdata[index[i]], ydata[index[i]], xdata[index[i+1]], ydata[index[i+1]], xdata[index[i]+j], ydata[index[i]+j]), 2); /*printf("dist = %g\n", dist);*/ } i++; } while(index[i] != n-1); return dist/(double)n; } double nsl_geom_linesim_area_error(const double xdata[], const double ydata[], const size_t n, const size_t index[]) { double area = 0; size_t i = 0, j; /* i: index of index[] */ do { /* for every point not in index[] calculate area */ /*printf("i=%d (index[i]-index[i+1]=%d-%d)\n", i , index[i], index[i+1]);*/ for (j = 1; j < index[i+1]-index[i]; j++) { /*printf("j=%d: area between %d %d %d\n", j, index[i]+j-1, index[i]+j, index[i+1]);*/ area += nsl_geom_three_point_area(xdata[index[i]+j-1], ydata[index[i]+j-1], xdata[index[i]+j], ydata[index[i]+j], xdata[index[i+1]], ydata[index[i+1]]); /*printf("area = %g\n", area);*/ } i++; } while(index[i] != n-1); return area/(double)n; } double nsl_geom_linesim_clip_diag_perpoint(const double xdata[], const double ydata[], const size_t n) { double dx = nsl_stats_maximum(xdata, n, NULL) - nsl_stats_minimum(xdata, n, NULL); double dy = nsl_stats_maximum(ydata, n, NULL) - nsl_stats_minimum(ydata, n, NULL); double d = sqrt(dx*dx+dy*dy); return d/(double)n; /* per point */ } double nsl_geom_linesim_clip_area_perpoint(const double xdata[], const double ydata[], const size_t n) { double dx = nsl_stats_maximum(xdata, n, NULL) - nsl_stats_minimum(xdata, n, NULL); double dy = nsl_stats_maximum(ydata, n, NULL) - nsl_stats_minimum(ydata, n, NULL); double A = dx*dy; return A/(double)n; /* per point */ } double nsl_geom_linesim_avg_dist_perpoint(const double xdata[], const double ydata[], const size_t n) { double dist = 0, dx, dy; size_t i; for (i = 0; i < n-1; i++) { dx = xdata[i+1]-xdata[i]; dy = ydata[i+1]-ydata[i]; dist += sqrt(dx*dx+dy*dy); } dist /= (double)n; return dist; } /*********** simplification algorithms *********/ void nsl_geom_linesim_douglas_peucker_step(const double xdata[], const double ydata[], const size_t start, const size_t end, size_t *nout, const double tol, size_t index[]) { /*printf("DP: %d - %d\n", start, end);*/ size_t i, nkey = start; double dist, maxdist = 0; /* search for key (biggest perp. distance) */ for (i = start+1; i < end; i++) { dist = nsl_geom_point_line_dist(xdata[start], ydata[start], xdata[end], ydata[end], xdata[i], ydata[i]); if (dist > maxdist) { maxdist=dist; nkey=i; } } /*printf("maxdist = %g @ i = %zu\n", maxdist, nkey);*/ if (maxdist > tol) { /*printf("take %d\n", nkey);*/ index[(*nout)++]=nkey; if (nkey-start > 1) nsl_geom_linesim_douglas_peucker_step(xdata, ydata, start, nkey, nout, tol, index); if (end-nkey > 1) nsl_geom_linesim_douglas_peucker_step(xdata, ydata, nkey, end, nout, tol, index); } /*printf("nout=%d\n", *nout);*/ } size_t nsl_geom_linesim_douglas_peucker(const double xdata[], const double ydata[], const size_t n, const double tol, size_t index[]) { size_t nout = 0; /*first point*/ index[nout++] = 0; nsl_geom_linesim_douglas_peucker_step(xdata, ydata, 0, n-1, &nout, tol, index); /* last point */ if (index[nout-1] != n-1) index[nout++] = n-1; /* sort array index */ nsl_sort_size_t(index, nout); return nout; } size_t nsl_geom_linesim_douglas_peucker_auto(const double xdata[], const double ydata[], const size_t n, size_t index[]) { double tol = nsl_geom_linesim_clip_diag_perpoint(xdata, ydata, n); return nsl_geom_linesim_douglas_peucker(xdata, ydata, n, tol, index); } double nsl_geom_linesim_douglas_peucker_variant(const double xdata[], const double ydata[], const size_t n, const size_t nout, size_t index[]) { size_t i; if (nout >= n) { /* all points */ for (i = 0; i < n; i++) index[i] = i; return 0; } /* first and last point */ size_t ntmp = 0; index[ntmp++] = 0; index[ntmp++] = n-1; if (nout <= 2) /* using first and last point */ return DBL_MAX; double *dist = (double *)malloc(n * sizeof(double)); double *maxdist = (double *)malloc(nout * sizeof(double)); /* max dist per edge */ + if (dist == NULL || maxdist == NULL) { + printf("nsl_geom_linesim_douglas_peucker_variant(): ERROR allocating memory!\n"); + return 0; + } for (i = 0; i < n; i++) { /* initialize dist */ dist[i] = nsl_geom_point_line_dist(xdata[0], ydata[0], xdata[n-1], ydata[n-1], xdata[i], ydata[i]); /*printf("%zu: %g\n", i, dist[i]);*/ } for (i = 0; i < nout; i++) maxdist[i] = 0; double newmaxdist; while (ntmp < nout) { size_t key = 0, v; /* find edge of maximum */ size_t maxindex; nsl_stats_maximum(maxdist, ntmp, &maxindex); /*printf("found edge of max at index %zu\n", maxindex);*/ /*newmaxdist = nsl_stats_maximum(dist, n, &key);*/ newmaxdist = 0; for (i = index[maxindex]+1; i < index[maxindex+1]; i++) { /*printf("i=%zu\n", i);*/ if (dist[i] > newmaxdist) { newmaxdist = dist[i]; key = i; } } /*printf("found key %zu (dist = %g)\n", key, newmaxdist);*/ ntmp++; dist[key] = 0; /* find index of previous key */ size_t previndex = 0; while (index[previndex+1] < key) previndex++; /*printf("previndex = %zu (update key %zu - %zu)\n", previndex, index[previndex], index[previndex+1]);*/ /* shift maxdist */ for (v = ntmp; v > previndex; v--) maxdist[v] = maxdist[v-1]; /* update dist[]. no update on last key */ if (ntmp < nout) { double tmpmax = 0; for (v = index[previndex]+1; v < key; v++) { /*printf("%zu to %zu - %zu", v, index[previndex], key);*/ dist[v] = nsl_geom_point_line_dist(xdata[index[previndex]], ydata[index[previndex]], xdata[key], ydata[key], xdata[v], ydata[v]); if (dist[v] > tmpmax) tmpmax = dist[v]; /*printf(" dist = %g\n", dist[v]);*/ } maxdist[previndex] = tmpmax; tmpmax = 0; for (v = key+1; v < index[previndex+1]; v++) { /*printf("%zu to %zu - %zu", v, key, index[previndex+1]);*/ dist[v] = nsl_geom_point_line_dist(xdata[key], ydata[key], xdata[index[previndex+1]], ydata[index[previndex+1]], xdata[v], ydata[v]); if (dist[v] > tmpmax) tmpmax = dist[v]; /*printf(" dist = %g\n", dist[v]);*/ } maxdist[previndex+1] = tmpmax; } /* put into index array */ for (v = ntmp; v > previndex+1; v--) index[v] = index[v-1]; index[previndex+1] = key; } free(dist); free(maxdist); return newmaxdist; } size_t nsl_geom_linesim_nthpoint(const size_t n, const size_t step, size_t index[]) { if (step < 1) { printf("step size must be > 0 (given: %zd)\n", step); return 0; } size_t i, nout = 0; /*first point*/ index[nout++] = 0; for (i = 1; i < n-1; i++) if (i%step == 0) index[nout++] = i; /* last point */ index[nout++] = n-1; return nout; } size_t nsl_geom_linesim_raddist(const double xdata[], const double ydata[], const size_t n, const double tol, size_t index[]) { size_t i, nout = 0, key = 0; /*first point*/ index[nout++] = 0; for (i = 1; i < n-1; i++) { /* distance to key point */ double dist = nsl_geom_point_point_dist(xdata[i], ydata[i], xdata[key], ydata[key]); /* distance to last point */ double lastdist = nsl_geom_point_point_dist(xdata[i], ydata[i], xdata[n-1], ydata[n-1]); /*printf("%d: %g %g\n", i, dist, lastdist);*/ if (dist > tol && lastdist > tol) { index[nout++] = i; key = i; } } /* last point */ index[nout++] = n-1; return nout; } size_t nsl_geom_linesim_raddist_auto(const double xdata[], const double ydata[], const size_t n, size_t index[]) { double tol = 10.*nsl_geom_linesim_clip_diag_perpoint(xdata, ydata, n); return nsl_geom_linesim_raddist(xdata, ydata, n, tol, index); } size_t nsl_geom_linesim_perpdist(const double xdata[], const double ydata[], const size_t n, const double tol, size_t index[]) { size_t nout = 0, key = 0, i; /*first point*/ index[nout++] = 0; for (i = 1; i < n-1; i++) { /* distance of point i to line key -- i+1 */ double dist = nsl_geom_point_line_dist(xdata[key], ydata[key], xdata[i+1], ydata[i+1], xdata[i], ydata[i]); /*printf("%d: %g\n", i, dist);*/ if (dist > tol) { /* take it */ index[nout++] = i; key = i; /*printf("%d: take it (key = %d)\n", i, key);*/ } else { /* ignore it */ if (i+1 < n-1) /* last point is taken anyway */ index[nout++] = i+1; /* take next point in any case */ /*printf("%d: ignore it (key = %d)\n", i, i+1);*/ key = ++i; } } /* last point */ index[nout++] = n-1; return nout; } size_t nsl_geom_linesim_perpdist_auto(const double xdata[], const double ydata[], const size_t n, size_t index[]) { double tol = nsl_geom_linesim_clip_diag_perpoint(xdata, ydata, n); return nsl_geom_linesim_perpdist(xdata, ydata, n, tol, index); } size_t nsl_geom_linesim_perpdist_repeat(const double xdata[], const double ydata[], const size_t n, const double tol, const size_t repeat, size_t index[]) { size_t i, j, nout; double *xtmp = (double *) malloc(n*sizeof(double)); double *ytmp = (double *) malloc(n*sizeof(double)); size_t *tmpindex = (size_t *) malloc(n*sizeof(size_t)); + if (xtmp == NULL || ytmp == NULL || tmpindex == NULL) { + printf("nsl_geom_linesim_perpdist_repeat(): ERROR allocating memory!\n"); + return 0; + } /* first round */ nout = nsl_geom_linesim_perpdist(xdata, ydata, n, tol, index); /* repeats */ for (i = 0; i < repeat - 1; i++) { for (j = 0; j < nout; j++) { xtmp[j] = xdata[index[j]]; ytmp[j] = ydata[index[j]]; tmpindex[j]= index[j]; /*printf("%g %g\n", xtmp[j], ytmp[j]);*/ } size_t tmpnout = nsl_geom_linesim_perpdist(xtmp, ytmp, nout, tol, tmpindex); for (j = 0; j < tmpnout; j++) { index[j] = index[tmpindex[j]]; /*printf("tmpindex[%d]: %d\n", j, tmpindex[j]);*/ } if (tmpnout == nout) /* return if nout does not change anymore */ break; else nout = tmpnout; } free(tmpindex); free(xtmp); free(ytmp); return nout; } size_t nsl_geom_linesim_interp(const double xdata[], const double ydata[], const size_t n, const double tol, size_t index[]) { size_t i, nout = 0; /*first point*/ index[nout++] = 0; size_t key = 0; for (i = 1; i < n-1; i++) { /*printf("%d: %d-%d\n", i, key, i+1);*/ double dist = nsl_geom_point_line_dist_y(xdata[key], ydata[key], xdata[i+1], ydata[i+1], xdata[i], ydata[i]); /*printf("%d: dist = %g\n", i, dist);*/ if (dist > tol) { /*printf("take it %d\n", i);*/ index[nout++] = key = i; } } /* last point */ index[nout++] = n-1; return nout; } size_t nsl_geom_linesim_interp_auto(const double xdata[], const double ydata[], const size_t n, size_t index[]) { double tol = nsl_geom_linesim_clip_diag_perpoint(xdata, ydata, n); return nsl_geom_linesim_interp(xdata, ydata, n, tol, index); } size_t nsl_geom_linesim_visvalingam_whyatt(const double xdata[], const double ydata[], const size_t n, const double tol, size_t index[]) { if (n < 3) /* we need at least three points */ return 0; size_t i, nout = n; double *area = (double *) malloc((n-2)*sizeof(double)); /* area associated with every point */ + if (area == NULL) { + printf("nsl_geom_linesim_visvalingam_whyatt(): ERROR allocating memory!\n"); + return 0; + } for (i = 0; i < n; i++) index[i] = i; for (i = 1; i < n-1; i++) area[i-1] = nsl_geom_three_point_area(xdata[i-1], ydata[i-1], xdata[i], ydata[i], xdata[i+1], ydata[i+1]); double minarea; size_t minindex; while ( (minarea = nsl_stats_minimum(area, n-2, &minindex)) < tol && nout > 2) { /*for (i=0; i < n-2; i++) if (area[i] 0) before--; while (index[after] == 0 && after < n-1) after++; if (minindex > 0) { /*before */ if (before > 0) { size_t beforebefore=before-1; while (index[beforebefore] == 0 && beforebefore > 0) beforebefore--; /*printf("recalculate area[%zu] from %zu %zu %zu\n", before-1, beforebefore, before, after);*/ tmparea = nsl_geom_three_point_area(xdata[beforebefore], ydata[beforebefore], xdata[before], ydata[before], xdata[after], ydata[after]); if (tmparea > area[before-1]) /* take largest value of new and old area */ area[before-1] = tmparea; } } if (minindex < n-3) { /* after */ if (after < n-1) { size_t afterafter = after+1; while (index[afterafter] == 0 && afterafter < n-1) afterafter++; /*printf("recalculate area[%zu] from %zu %zu %zu\n",after-1, before, after, afterafter);*/ tmparea = nsl_geom_three_point_area(xdata[before], ydata[before], xdata[after], ydata[after], xdata[afterafter], ydata[afterafter]); if (tmparea > area[after-1]) /* take largest value of new and old area */ area[after-1] = tmparea; } } nout--; }; /*for(i=0;i tol) { /* take it */ /*printf("%d: take it\n", i);*/ key = i-1; key2 = i; index[nout++] = i-1; } } /* last point */ index[nout++] = n-1; return nout; } size_t nsl_geom_linesim_reumann_witkam_auto(const double xdata[], const double ydata[], const size_t n, size_t index[]) { double tol = nsl_geom_linesim_clip_diag_perpoint(xdata, ydata, n); return nsl_geom_linesim_reumann_witkam(xdata, ydata, n, tol, index); } size_t nsl_geom_linesim_opheim(const double xdata[], const double ydata[], const size_t n, const double mintol, const double maxtol, size_t index[]) { size_t i, nout = 0, key = 0, key2; /*first point*/ index[nout++] = 0; for (i = 1; i < n-1; i++) { double dist, perpdist; do { /* find key2 */ dist = nsl_geom_point_point_dist(xdata[key], ydata[key], xdata[i], ydata[i]); /*printf("dist = %g: %d-%d\n", dist, key, i);*/ i++; } while (dist < mintol); i--; if (key == i-1) /*i+1 outside mintol */ key2 = i; else key2 = i-1; /* last point inside */ /*printf("found key2 @%d\n", key2);*/ do { /* find next key */ dist = nsl_geom_point_point_dist(xdata[key], ydata[key], xdata[i], ydata[i]); perpdist = nsl_geom_point_line_dist(xdata[key], ydata[key], xdata[key2], ydata[key2], xdata[i], ydata[i]); /*printf("dist = %g, perpdist=%g: %d\n", dist, perpdist, i);*/ i++; } while (dist < maxtol && perpdist < mintol); i--; if (key == i-1) /*i+1 outside */ key = i; else { key = i-1; i--; } /*printf("new key: %d\n", key);*/ index[nout++] = key; } /* last point */ if (index[nout-1] != n-1) index[nout++] = n-1; return nout; } size_t nsl_geom_linesim_opheim_auto(const double xdata[], const double ydata[], const size_t n, size_t index[]) { double mintol = 10.*nsl_geom_linesim_clip_diag_perpoint(xdata, ydata, n); /* TODO: calculate max tolerance ? */ double maxtol = 5.*mintol; return nsl_geom_linesim_opheim(xdata, ydata, n, mintol, maxtol, index); return 0; } size_t nsl_geom_linesim_lang(const double xdata[], const double ydata[], const size_t n, const double tol, const size_t region, size_t index[]) { size_t i, j, nout = 0, key = 0; /*first point*/ index[nout++] = 0; double dist, maxdist; for (i = 1; i < n-1; i++) { size_t tmpregion=region; if (key+tmpregion > n-1) /* end of data set */ tmpregion = n-1-key; do { maxdist = 0; for (j = 1; j < tmpregion; j++) { dist = nsl_geom_point_line_dist(xdata[key], ydata[key], xdata[key+tmpregion], ydata[key+tmpregion], xdata[key+j], ydata[key+j]); /*printf("%d: dist (%d to %d-%d) = %g\n", j, key+j, key, key+tmpregion, dist);*/ if (dist > maxdist) maxdist = dist; } /*printf("tol = %g maxdist = %g\n", tol, maxdist);*/ tmpregion--; /*printf("region = %d\n", tmpregion);*/ } while (maxdist>tol && tmpregion>0); i += tmpregion; index[nout++] = key = i; /*printf("take it (%d) key=%d\n", i, key);*/ } /* last point */ if (index[nout-1] != n-1) index[nout++] = n-1; return nout; } size_t nsl_geom_linesim_lang_auto(const double xdata[], const double ydata[], const size_t n, size_t index[]) { double tol = nsl_geom_linesim_clip_diag_perpoint(xdata, ydata, n); /* TODO: calculate search region */ double region = 0; printf("nsl_geom_linesim_lang_auto(): Not implemented yet\n"); return nsl_geom_linesim_lang(xdata, ydata, n, tol, region, index); }