diff --git a/dspl/src/gnuplot/contour2d.c b/dspl/src/gnuplot/contour2d.c
new file mode 100644
index 0000000..4bc2393
--- /dev/null
+++ b/dspl/src/gnuplot/contour2d.c
@@ -0,0 +1,558 @@
+/*
+* Copyright (c) 2015-2024 Sergey Bakhurin
+* Digital Signal Processing Library [http://dsplib.org]
+*
+* This file is part of DSPL.
+*
+* 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.
+*
+* DSPL 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 Foobar. If not, see .
+*/
+#include
+#include
+#include
+#include "dspl.h"
+
+
+
+#define BSIZE 256
+#define DIST_TH 1E-8
+
+typedef struct
+{
+ point2d_t p[2];
+ int flag;
+} linseg_t;
+
+
+
+int add_linseg(linseg_t** ls, int* lsnum, int* lscnt,
+ point2d_t* p0, point2d_t* p1)
+
+{
+ int n, c;
+ n = *lsnum;
+ c = *lscnt;
+ // проверяем выделение памяти.
+ if((n == 0) && ((*ls)==NULL))
+ {
+ n = BSIZE;
+ (*ls) = (linseg_t*)malloc(n * sizeof(linseg_t));
+ }
+ else
+ {
+ // при необходимости увеличиваем
+ if(c >= n)
+ {
+ n += BSIZE;
+ (*ls) = (linseg_t*)realloc((*ls), n * sizeof(linseg_t));
+ }
+ }
+ (*ls)[c].p[0][0] = p0[0][0];
+ (*ls)[c].p[0][1] = p0[0][1];
+ (*ls)[c].p[1][0] = p1[0][0];
+ (*ls)[c].p[1][1] = p1[0][1];
+ (*ls)[c].flag = 1;
+ c++;
+ (*lsnum) = n;
+ (*lscnt) = c;
+
+ return RES_OK;
+}
+
+int linseg_create(double* z, double* x, double* y,
+ int n, int m, double lev,
+ linseg_t** ls, int* sz)
+{
+ int lsnum, lscnt, t, in, im, i;
+ point2d_t p0 = {0};
+ point2d_t p1 = {0};
+
+ double dx;
+ double dy;
+
+ if((ls== NULL)||(z==NULL))
+ return ERROR_PTR;
+ lsnum = 0;
+ lscnt = 0;
+
+ for(in = 0; in < n-1; in++)
+ {
+ for(im = 0; im < m-1; im++)
+ {
+ i = in + im * n;
+ t = 0;
+ t += z[i] > lev ? 8 : 0;
+ t += z[i+n] > lev ? 4 : 0;
+ t += z[i+n+1] > lev ? 2 : 0;
+ t += z[i+1] > lev ? 1 : 0;
+
+ printf("%d, %d, %d\n", in, im, t);
+ switch(t)
+ {
+ case 0:
+ case 15:
+ break;
+ case 1:
+ case 14:
+ dx = (lev - z[i]) / (z[i+1] - z[i]);
+ p0[0] = x[in] + dx;
+ p0[1] = y[im];
+
+ dy = (lev - z[i+1]) / (z[i+n+1] - z[i+1]);
+ p1[0] = x[in+1];
+ p1[1] = y[im] + dy;
+ add_linseg(ls, &lsnum, &lscnt, &p0, &p1);
+ break;
+ case 2:
+ case 13:
+ dx = (lev - z[i+n]) / (z[i+n+1] - z[i+n]);
+ p0[0] = x[in] + dx;
+ p0[1] = y[im + 1];
+
+ dy = (lev - z[i+1]) / (z[i+n+1] - z[i+1]);
+ p1[0] = x[in + 1];
+ p1[1] = y[im] + dy;
+ add_linseg(ls, &lsnum, &lscnt, &p0, &p1);
+ break;
+ case 3:
+ case 12:
+ dx = (lev - z[i]) / (z[i+1] - z[i]);
+ p0[0] = x[in] + dx;
+ p0[1] = y[im];
+
+ dx = (lev - z[i+n]) / (z[i+n+1] - z[i+n]);
+ p1[0] = x[in] + dx;
+ p1[1] = y[im+1];
+ add_linseg(ls, &lsnum, &lscnt, &p0, &p1);
+ break;
+ case 4:
+ case 11:
+ dy = (lev - z[i]) / (z[i+n] - z[i]);
+ p0[0] = x[in];
+ p0[1] = y[im] + dy;
+
+ dx = (lev - z[i+n]) / (z[i+n+1] - z[i+n]);
+ p1[0] = x[in] + dx;
+ p1[1] = y[im + 1];
+ add_linseg(ls, &lsnum, &lscnt, &p0, &p1);
+ break;
+ case 5:
+ dy = (lev - z[i]) / (z[i+n] - z[i]);
+ p0[0] = x[in];
+ p0[1] = y[im] + dy;
+ dx = (lev - z[i]) / (z[i+1] - z[i]);
+ p1[0] = x[in] + dx;
+ p1[1] = y[im];
+ add_linseg(ls, &lsnum, &lscnt, &p0, &p1);
+
+ dx = (lev - z[i+n]) / (z[i+n+1] - z[i+n]);
+ p0[0] = x[in]+dx;
+ p0[1] = y[im+1];
+ dy = (lev - z[i+1]) / (z[i+n+1] - z[i+1]);
+ p1[0] = x[in+1];
+ p1[1] = y[im] + dy;
+ add_linseg(ls, &lsnum, &lscnt, &p0, &p1);
+ break;
+ case 6:
+ case 9:
+ dy = (lev - z[i]) / (z[i+n] - z[i]);
+ p0[0] = x[in];
+ p0[1] = y[im] + dy;
+ dy = (lev - z[i+1]) / (z[i+n+1] - z[i+1]);
+ p1[0] = x[in+1];
+ p1[1] = y[im]+dy;
+ add_linseg(ls, &lsnum, &lscnt, &p0, &p1);
+ break;
+
+ case 7:
+ case 8:
+ dy = (lev - z[i]) / (z[i+n] - z[i]);
+ p0[0] = x[in];
+ p0[1] = y[im] + dy;
+ dx = (lev - z[i]) / (z[i+1] - z[i]);
+ p1[0] = x[in] + dx;
+ p1[1] = y[im];
+ add_linseg(ls, &lsnum, &lscnt, &p0, &p1);
+ break;
+ case 10:
+ dy = (lev - z[i]) / (z[i+n] - z[i]);
+ p0[0] = x[in];
+ p0[1] = y[im] + dy;
+ dx = (lev - z[i+n]) / (z[i+n+1] - z[i+1]);
+ p1[0] = x[in] + dx;
+ p1[1] = y[im+1];
+ add_linseg(ls, &lsnum, &lscnt, &p0, &p1);
+
+ dx = (lev - z[i]) / (z[i+1] - z[i]);
+ p0[0] = x[in]+dx;
+ p0[1] = y[im];
+ dy = (lev - z[i+1]) / (z[i+n+1] - z[i+1]);
+ p1[0] = x[in+1];
+ p1[1] = y[im] + dy;
+ add_linseg(ls, &lsnum, &lscnt, &p0, &p1);
+ break;
+ default:
+ break;
+ }
+ }
+ }
+ *ls = (linseg_t*)realloc(*ls, lscnt * sizeof(linseg_t));
+ *sz = lscnt;
+ return RES_OK;
+}
+
+double dist(point2d_t* p0, point2d_t* p1)
+{
+ double dx, dy;
+ dx = p0[0][0] - p1[0][0];
+ dy = p0[0][1] - p1[0][1];
+ return sqrt(dx*dx + dy*dy);
+}
+
+int line_create(linseg_t* ls, int nls, line2d_t* line)
+{
+ int i, j, c, n, k;
+ if(!line || !ls)
+ return ERROR_PTR;
+ // printf("*line = %x, *np = %d \n", *line, *np);
+ i = 0;
+ while(!(ls[i].flag) && i < nls)
+ i++;
+ if(i==nls)
+ {
+ line->npoints = 0;
+ return RES_OK;
+ }
+ //printf("i = %d, ls[%d] = [%.1f %.1f] -- [%.1f, %.1f]\n",
+ // i, i, ls[i].p[0][0], ls[i].p[0][1],
+ // ls[i].p[1][0], ls[i].p[1][1]);
+ n = BSIZE;
+ if(line->points == NULL)
+ line->points = (point2d_t*)malloc(n*sizeof(point2d_t));
+ else
+ line->points = (point2d_t*)realloc(line->points, n*sizeof(point2d_t));
+
+ c = 0;
+ line->points[c][0] = ls[i].p[0][0];
+ line->points[c][1] = ls[i].p[0][1];
+ c++;
+ line->points[c][0] = ls[i].p[1][0];
+ line->points[c][1] = ls[i].p[1][1];
+ c++;
+ ls[i].flag = 0;
+
+ for(i = 0; i < nls; i++)
+ {
+ for(j = 0; j < nls; j++)
+ {
+ if(ls[j].flag)
+ {
+ //сравниваем последнюю точку линии с первой точкой отрезка ls[j]
+ if(dist(line->points+c-1, ls[j].p) < DIST_TH)
+ {
+
+ // проверяем выделение памяти.
+ // при необходимости увеличиваем
+ if(c>=n)
+ {
+ n += BSIZE;
+ line->points = (point2d_t*)realloc(line->points, n*sizeof(point2d_t));
+ }
+ //если первая точка отрезка совпадает, то добавляем
+ //в линию вторую точку СЗАДИ
+ line->points[c][0] = ls[j].p[1][0];
+ line->points[c][1] = ls[j].p[1][1];
+ ls[j].flag = 0;
+ c++;
+ }
+ }
+ if(ls[j].flag)
+ {
+ // сравниваем последнюю точку линии
+ // со второй точкой отрезка ls[j]
+ if(dist(line->points+c-1, ls[j].p+1) < DIST_TH)
+ {
+ // проверяем выделение памяти.
+ // при необходимости увеличиваем
+ if(c>=n)
+ {
+ n += BSIZE;
+ line->points = (point2d_t*)realloc(line->points, n*sizeof(point2d_t));
+ }
+ //если вторая точка совпадает, то добавляем
+ //в линию первую точку СЗАДИ
+ line->points[c][0] = ls[j].p[0][0];
+ line->points[c][1] = ls[j].p[0][1];
+ ls[j].flag = 0;
+ c++;
+ }
+ }
+
+ if(ls[j].flag)
+ {
+ //сравниваем последнюю точку линии с первой точкой отрезка ls[j]
+ if(dist(line->points, ls[j].p) < DIST_TH)
+ {
+
+ // проверяем выделение памяти.
+ // при необходимости увеличиваем
+ if(c>=n)
+ {
+ n += BSIZE;
+ line->points = (point2d_t*)realloc(line->points, n*sizeof(point2d_t));
+ }
+ //если первая точка отрезка совпадает, то добавляем
+ //в линию вторую точку СПЕРЕДИ
+ for(k = c; k>0; k--)
+ {
+ line->points[k][0] = line->points[k-1][0];
+ line->points[k][1] = line->points[k-1][1];
+ }
+ line->points[0][0] = ls[j].p[1][0];
+ line->points[0][1] = ls[j].p[1][1];
+ ls[j].flag = 0;
+ c++;
+ }
+ }
+ if(ls[j].flag)
+ {
+ // сравниваем последнюю точку линии
+ // со второй точкой отрезка ls[j]
+ if(dist(line->points, ls[j].p+1) < DIST_TH)
+ {
+ // проверяем выделение памяти.
+ // при необходимости увеличиваем
+ if(c>=n)
+ {
+ n += BSIZE;
+ line->points = (point2d_t*)realloc(line->points, n*sizeof(point2d_t));
+ }
+ //если вторая точка совпадает, то добавляем
+ //в линию первую точку СПЕРЕДИ
+
+ for(k = c; k>0; k--)
+ {
+ line->points[k][0] = line->points[k-1][0];
+ line->points[k][1] = line->points[k-1][1];
+ }
+
+ line->points[0][0] = ls[j].p[0][0];
+ line->points[0][1] = ls[j].p[0][1];
+ ls[j].flag = 0;
+ c++;
+ }
+ }
+
+ }
+ }
+ line->points = (point2d_t*)realloc(line->points, c*sizeof(point2d_t));
+ line->npoints = c;
+ return RES_OK;
+}
+
+
+
+
+int DSPL_API contour2d(double* z, double* x, double* y,
+ int n, int m, double lev,
+ contour2d_t* c)
+{
+
+ linseg_t *ls = NULL;
+ int nls, i, j, nlines;
+
+
+
+ linseg_create(z, x, y, n, m, lev, &ls, &nls);
+ for(i =0; i< nls; i++)
+ printf("%d, [%.1f %.1f] -- [%.1f %.1f]\n", i, ls[i].p[0][0], ls[i].p[0][1], ls[i].p[1][0], ls[i].p[1][1]);
+
+ if (c->lines == NULL)
+ {
+ nlines = BSIZE;
+ c->lines = (line2d_t*)malloc(nlines*sizeof(line2d_t));
+ }
+
+
+ i = -1;
+ do
+ {
+ i++;
+ if(i >= nlines)
+ {
+ nlines+= BSIZE;
+ c->lines = (line2d_t*)realloc(c->lines, nlines*sizeof(line2d_t));
+ }
+ line_create(ls, nls, c->lines+i);
+ for(j =0; j< c->lines[i].npoints; j++)
+ printf("[%.1f %.1f] -- ", c->lines[i].points[j][0], c->lines[i].points[j][1]);
+ printf("\n");
+ }while(c->lines[i].npoints);
+ c->lines = (line2d_t*)realloc(c->lines, i*sizeof(line2d_t));
+ c->nlines = i;
+ c->level = lev;
+
+ if(ls)
+ free(ls);
+ return RES_OK;
+}
+
+
+
+
+/*
+int line_create(linseg_t* ls, int nls, point2d_t** line, int* np)
+{
+ int i, j, c, n, k;
+ if(!line || !ls || !np)
+ return ERROR_PTR;
+ // printf("*line = %x, *np = %d \n", *line, *np);
+ i = 0;
+ while(!(ls[i].flag) && i < nls)
+ i++;
+ if(i==nls)
+ {
+ *np = 0;
+ return RES_OK;
+ }
+ //printf("i = %d, ls[%d] = [%.1f %.1f] -- [%.1f, %.1f]\n",
+ // i, i, ls[i].p[0][0], ls[i].p[0][1],
+ // ls[i].p[1][0], ls[i].p[1][1]);
+ n = BSIZE;
+ if((*line) == NULL)
+ (*line) = (point2d_t*)malloc(n*sizeof(point2d_t));
+ else
+ (*line) = (point2d_t*)realloc((*line), n*sizeof(point2d_t));
+
+ c = 0;
+ (*line)[c][0] = ls[i].p[0][0];
+ (*line)[c][1] = ls[i].p[0][1];
+ c++;
+ (*line)[c][0] = ls[i].p[1][0];
+ (*line)[c][1] = ls[i].p[1][1];
+ c++;
+ ls[i].flag = 0;
+
+ for(i = 0; i < nls; i++)
+ {
+ for(j = 0; j < nls; j++)
+ {
+ if(ls[j].flag)
+ {
+ //сравниваем последнюю точку линии с первой точкой отрезка ls[j]
+ if(dist((*line)+c-1, ls[j].p) < DIST_TH)
+ {
+
+ // проверяем выделение памяти.
+ // при необходимости увеличиваем
+ if(c>=n)
+ {
+ n += BSIZE;
+ (*line) = (point2d_t*)realloc(*line, n*sizeof(point2d_t));
+ }
+ //если первая точка отрезка совпадает, то добавляем
+ //в линию вторую точку СЗАДИ
+ (*line)[c][0] = ls[j].p[1][0];
+ (*line)[c][1] = ls[j].p[1][1];
+ ls[j].flag = 0;
+ c++;
+ }
+ }
+ if(ls[j].flag)
+ {
+ // сравниваем последнюю точку линии
+ // со второй точкой отрезка ls[j]
+ if(dist((*line)+c-1, ls[j].p+1) < DIST_TH)
+ {
+ // проверяем выделение памяти.
+ // при необходимости увеличиваем
+ if(c>=n)
+ {
+ n += BSIZE;
+ (*line) = (point2d_t*)realloc(*line, n*sizeof(point2d_t));
+ }
+ //если вторая точка совпадает, то добавляем
+ //в линию первую точку СЗАДИ
+ (*line)[c][0] = ls[j].p[0][0];
+ (*line)[c][1] = ls[j].p[0][1];
+ ls[j].flag = 0;
+ c++;
+ }
+ }
+
+ if(ls[j].flag)
+ {
+ //сравниваем последнюю точку линии с первой точкой отрезка ls[j]
+ if(dist((*line), ls[j].p) < DIST_TH)
+ {
+
+ // проверяем выделение памяти.
+ // при необходимости увеличиваем
+ if(c>=n)
+ {
+ n += BSIZE;
+ (*line) = (point2d_t*)realloc(*line, n*sizeof(point2d_t));
+ }
+ //если первая точка отрезка совпадает, то добавляем
+ //в линию вторую точку СПЕРЕДИ
+ for(k = c; k>0; k--)
+ {
+ (*line)[k][0] = (*line)[k-1][0];
+ (*line)[k][1] = (*line)[k-1][1];
+ }
+ (*line)[0][0] = ls[j].p[1][0];
+ (*line)[0][1] = ls[j].p[1][1];
+ ls[j].flag = 0;
+ c++;
+ }
+ }
+ if(ls[j].flag)
+ {
+ // сравниваем последнюю точку линии
+ // со второй точкой отрезка ls[j]
+ if(dist((*line), ls[j].p+1) < DIST_TH)
+ {
+ // проверяем выделение памяти.
+ // при необходимости увеличиваем
+ if(c>=n)
+ {
+ n += BSIZE;
+ (*line) = (point2d_t*)realloc(*line, n*sizeof(point2d_t));
+ }
+ //если вторая точка совпадает, то добавляем
+ //в линию первую точку СПЕРЕДИ
+
+ for(k = c; k>0; k--)
+ {
+ (*line)[k][0] = (*line)[k-1][0];
+ (*line)[k][1] = (*line)[k-1][1];
+ }
+
+ (*line)[0][0] = ls[j].p[0][0];
+ (*line)[0][1] = ls[j].p[0][1];
+ ls[j].flag = 0;
+ c++;
+ }
+ }
+
+ }
+ }
+ (*line) = (point2d_t*)realloc((*line), c*sizeof(point2d_t));
+ *np = c;
+ return RES_OK;
+}
+
+*/
+
+
+
diff --git a/examples/src/contour_test.c b/examples/src/contour_test.c
index 70da6f5..89d5127 100644
--- a/examples/src/contour_test.c
+++ b/examples/src/contour_test.c
@@ -5,351 +5,6 @@
#define N 3
#define M 4
-#define BSIZE 256
-
-typedef double point2d_t[2];
-
-typedef struct
-{
- point2d_t p[2];
- int flag;
-} linseg_t;
-
-
-int add_linseg(linseg_t** ls, int* lsnum, int* lscnt,
- point2d_t* p0, point2d_t* p1)
-
-{
- int n, c;
- n = *lsnum;
- c = *lscnt;
- // проверяем выделение памяти.
- if((n == 0) && ((*ls)==NULL))
- {
- n = BSIZE;
- (*ls) = (linseg_t*)malloc(n * sizeof(linseg_t));
- }
- else
- {
- // при необходимости увеличиваем
- if(c >= n)
- {
- n += BSIZE;
- (*ls) = (linseg_t*)realloc((*ls), n * sizeof(linseg_t));
- }
- }
- (*ls)[c].p[0][0] = p0[0][0];
- (*ls)[c].p[0][1] = p0[0][1];
- (*ls)[c].p[1][0] = p1[0][0];
- (*ls)[c].p[1][1] = p1[0][1];
- (*ls)[c].flag = 1;
- c++;
- (*lsnum) = n;
- (*lscnt) = c;
-
- return RES_OK;
-}
-
-int linseg_create(double* z, double* x, double* y,
- int n, int m, double lev,
- linseg_t** ls, int* sz)
-{
- int lsnum, lscnt, t, in, im, i;
- point2d_t p0 = {0};
- point2d_t p1 = {0};
-
- double dx;
- double dy;
-
- if((ls== NULL)||(z==NULL))
- return ERROR_PTR;
- lsnum = 0;
- lscnt = 0;
-
- for(in = 0; in < n-1; in++)
- {
- for(im = 0; im < m-1; im++)
- {
- i = in + im * n;
- t = 0;
- t += z[i] > lev ? 8 : 0;
- t += z[i+n] > lev ? 4 : 0;
- t += z[i+n+1] > lev ? 2 : 0;
- t += z[i+1] > lev ? 1 : 0;
-
- printf("%d, %d, %d\n", in, im, t);
- switch(t)
- {
- case 0:
- case 15:
- break;
- case 1:
- case 14:
- // z[i,j] * (1-dx) + z[i+1, j] * dx = lev
- // z[i,j] + dx *(z[i+1, j] - z[i,j]) = lev
- dx = (lev - z[i]) / (z[i+1] - z[i]);
- p0[0] = x[in] + dx;
- p0[1] = y[im];
-
- dy = (lev - z[i+1]) / (z[i+n+1] - z[i+1]);
- p1[0] = x[in+1];
- p1[1] = y[im] + dy;
- add_linseg(ls, &lsnum, &lscnt, &p0, &p1);
- break;
- case 2:
- case 13:
- dx = (lev - z[i+n]) / (z[i+n+1] - z[i+n]);
- p0[0] = x[in] + dx;
- p0[1] = y[im + 1];
-
- dy = (lev - z[i+1]) / (z[i+n+1] - z[i+1]);
- p1[0] = x[in + 1];
- p1[1] = y[im] + dy;
- add_linseg(ls, &lsnum, &lscnt, &p0, &p1);
- break;
- case 3:
- case 12:
- dx = (lev - z[i]) / (z[i+1] - z[i]);
- p0[0] = x[in] + dx;
- p0[1] = y[im];
-
- dx = (lev - z[i+n]) / (z[i+n+1] - z[i+n]);
- p1[0] = x[in] + dx;
- p1[1] = y[im+1];
- add_linseg(ls, &lsnum, &lscnt, &p0, &p1);
- break;
- case 4:
- case 11:
- dy = (lev - z[i]) / (z[i+n] - z[i]);
- p0[0] = x[in];
- p0[1] = y[im] + dy;
-
- dx = (lev - z[i+n]) / (z[i+n+1] - z[i+n]);
- p1[0] = x[in] + dx;
- p1[1] = y[im + 1];
- add_linseg(ls, &lsnum, &lscnt, &p0, &p1);
- break;
- case 5:
- dy = (lev - z[i]) / (z[i+n] - z[i]);
- p0[0] = x[in];
- p0[1] = y[im] + dy;
- dx = (lev - z[i]) / (z[i+1] - z[i]);
- p1[0] = x[in] + dx;
- p1[1] = y[im];
- add_linseg(ls, &lsnum, &lscnt, &p0, &p1);
-
- dx = (lev - z[i+n]) / (z[i+n+1] - z[i+n]);
- p0[0] = x[in]+dx;
- p0[1] = y[im+1];
- dy = (lev - z[i+1]) / (z[i+n+1] - z[i+1]);
- p1[0] = x[in+1];
- p1[1] = y[im] + dy;
- add_linseg(ls, &lsnum, &lscnt, &p0, &p1);
- break;
- case 6:
- case 9:
- dy = (lev - z[i]) / (z[i+n] - z[i]);
- p0[0] = x[in];
- p0[1] = y[im] + dy;
- dy = (lev - z[i+1]) / (z[i+n+1] - z[i+1]);
- p1[0] = x[in+1];
- p1[1] = y[im]+dy;
- add_linseg(ls, &lsnum, &lscnt, &p0, &p1);
- break;
-
- case 7:
- case 8:
- dy = (lev - z[i]) / (z[i+n] - z[i]);
- p0[0] = x[in];
- p0[1] = y[im] + dy;
- dx = (lev - z[i]) / (z[i+1] - z[i]);
- p1[0] = x[in] + dx;
- p1[1] = y[im];
- add_linseg(ls, &lsnum, &lscnt, &p0, &p1);
- break;
- case 10:
- dy = (lev - z[i]) / (z[i+n] - z[i]);
- p0[0] = x[in];
- p0[1] = y[im] + dy;
- dx = (lev - z[i+n]) / (z[i+n+1] - z[i+1]);
- p1[0] = x[in] + dx;
- p1[1] = y[im+1];
- add_linseg(ls, &lsnum, &lscnt, &p0, &p1);
-
- dx = (lev - z[i]) / (z[i+1] - z[i]);
- p0[0] = x[in]+dx;
- p0[1] = y[im];
- dy = (lev - z[i+1]) / (z[i+n+1] - z[i+1]);
- p1[0] = x[in+1];
- p1[1] = y[im] + dy;
- add_linseg(ls, &lsnum, &lscnt, &p0, &p1);
- break;
- default:
- break;
- }
- }
- }
- *ls = (linseg_t*)realloc(*ls, lscnt * sizeof(linseg_t));
- *sz = lscnt;
- return RES_OK;
-}
-
-
-double dist(point2d_t* p0, point2d_t* p1)
-{
- double dx, dy;
- dx = p0[0][0] - p1[0][0];
- dy = p0[0][1] - p1[0][1];
- return sqrt(dx*dx + dy*dy);
-}
-
-
-
-int line_create(linseg_t* ls, int nls, point2d_t** line, int* np)
-{
- int i, j, c, n, k;
- if(!line || !ls || !np)
- return ERROR_PTR;
- // printf("*line = %x, *np = %d \n", *line, *np);
- i = 0;
- while(!(ls[i].flag) && i < nls)
- i++;
- if(i==nls)
- {
- *np = 0;
- return RES_OK;
- }
- //printf("i = %d, ls[%d] = [%.1f %.1f] -- [%.1f, %.1f]\n",
- // i, i, ls[i].p[0][0], ls[i].p[0][1],
- // ls[i].p[1][0], ls[i].p[1][1]);
- n = BSIZE;
- if((*line) == NULL)
- (*line) = (point2d_t*)malloc(n*sizeof(point2d_t));
- else
- (*line) = (point2d_t*)realloc((*line), n*sizeof(point2d_t));
-
- c = 0;
- (*line)[c][0] = ls[i].p[0][0];
- (*line)[c][1] = ls[i].p[0][1];
- c++;
- (*line)[c][0] = ls[i].p[1][0];
- (*line)[c][1] = ls[i].p[1][1];
- c++;
- ls[i].flag = 0;
-
- for(i = 0; i < nls; i++)
- {
- for(j = 0; j < nls; j++)
- {
- if(ls[j].flag)
- {
- //сравниваем последнюю точку линии с первой точкой отрезка ls[j]
- if(dist((*line)+c-1, ls[j].p) < 1E-8)
- {
-
- // проверяем выделение памяти.
- // при необходимости увеличиваем
- if(c>=n)
- {
- n += BSIZE;
- (*line) = (point2d_t*)realloc(*line, n*sizeof(point2d_t));
- }
- //если первая точка отрезка совпадает, то добавляем
- //в линию вторую точку СЗАДИ
- (*line)[c][0] = ls[j].p[1][0];
- (*line)[c][1] = ls[j].p[1][1];
- ls[j].flag = 0;
- c++;
- }
- }
- if(ls[j].flag)
- {
- // сравниваем последнюю точку линии
- // со второй точкой отрезка ls[j]
- if(dist((*line)+c-1, ls[j].p+1) < 1E-8)
- {
- // проверяем выделение памяти.
- // при необходимости увеличиваем
- if(c>=n)
- {
- n += BSIZE;
- (*line) = (point2d_t*)realloc(*line, n*sizeof(point2d_t));
- }
- //если вторая точка совпадает, то добавляем
- //в линию первую точку СЗАДИ
- (*line)[c][0] = ls[j].p[0][0];
- (*line)[c][1] = ls[j].p[0][1];
- ls[j].flag = 0;
- c++;
- }
- }
-
- if(ls[j].flag)
- {
- //сравниваем последнюю точку линии с первой точкой отрезка ls[j]
- if(dist((*line), ls[j].p) < 1E-8)
- {
-
- // проверяем выделение памяти.
- // при необходимости увеличиваем
- if(c>=n)
- {
- n += BSIZE;
- (*line) = (point2d_t*)realloc(*line, n*sizeof(point2d_t));
- }
- //если первая точка отрезка совпадает, то добавляем
- //в линию вторую точку СПЕРЕДИ
- for(k = c; k>0; k--)
- {
- (*line)[k][0] = (*line)[k-1][0];
- (*line)[k][1] = (*line)[k-1][1];
- }
- (*line)[0][0] = ls[j].p[1][0];
- (*line)[0][1] = ls[j].p[1][1];
- ls[j].flag = 0;
- c++;
- }
- }
- if(ls[j].flag)
- {
- // сравниваем последнюю точку линии
- // со второй точкой отрезка ls[j]
- if(dist((*line), ls[j].p+1) < 1E-8)
- {
- // проверяем выделение памяти.
- // при необходимости увеличиваем
- if(c>=n)
- {
- n += BSIZE;
- (*line) = (point2d_t*)realloc(*line, n*sizeof(point2d_t));
- }
- //если вторая точка совпадает, то добавляем
- //в линию первую точку СПЕРЕДИ
-
- for(k = c; k>0; k--)
- {
- (*line)[k][0] = (*line)[k-1][0];
- (*line)[k][1] = (*line)[k-1][1];
- }
-
- (*line)[0][0] = ls[j].p[0][0];
- (*line)[0][1] = ls[j].p[0][1];
- ls[j].flag = 0;
- c++;
- }
- }
-
- }
- }
- (*line) = (point2d_t*)realloc((*line), c*sizeof(point2d_t));
- *np = c;
-
- return RES_OK;
-}
-
-
-
-
int main(int argc, char* argv[])
{
@@ -362,32 +17,13 @@ int main(int argc, char* argv[])
double z[N*M] = { 1, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1};
double x[N] = {0.0, 1.0, 2.0};
double y[M] = {0.0, 1.0, 2.0, 3.0};
- linseg_t *ls = NULL;
- int nls, i, j;
- point2d_t** line = NULL;
- int *np = NULL;
-
- line = (point2d_t**)malloc(10*sizeof(point2d_t*));
- memset(line,0,10*sizeof(point2d_t*));
-
- np = (int*)malloc(10*sizeof(int));
- memset(np, 0, 10*sizeof(int));
+ void* handle; /* DSPL handle */
+ handle = dspl_load(); /* Load DSPL function */
- linseg_create(z, x, y, N, M, 0.5, &ls, &nls);
- for(i =0; i< nls; i++)
- printf("%d, [%.1f %.1f] -- [%.1f %.1f]\n", i, ls[i].p[0][0], ls[i].p[0][1], ls[i].p[1][0], ls[i].p[1][1]);
-
-
- i = -1;
- do
- {
- i++;
- line_create(ls, nls, line+i, np+i);
- for(j =0; j< np[i]; j++)
- printf("[%.1f %.1f] -- ", line[i][j][0], line[i][j][1]);
- printf("\n");
- }while(np[i]);
+ contour2d_t contour = {0};
+ contour2d(z, x, y, N, M, 0.5,&contour);
+ dspl_free(handle); /* free dspl handle */
return 0;
}
diff --git a/include/dspl.c b/include/dspl.c
index 6daddf3..f9d82de 100644
--- a/include/dspl.c
+++ b/include/dspl.c
@@ -55,6 +55,7 @@ p_cheby2_ap_wp1 cheby2_ap_wp1 ;
p_cheby2_ap_zp cheby2_ap_zp ;
p_cmplx2re cmplx2re ;
p_concat concat ;
+p_contour2d contour2d ;
p_conv conv ;
p_conv_cmplx conv_cmplx ;
p_conv_fft conv_fft ;
@@ -274,6 +275,7 @@ void* dspl_load()
LOAD_FUNC(cheby2_ap_zp);
LOAD_FUNC(cmplx2re);
LOAD_FUNC(concat);
+ LOAD_FUNC(contour2d);
LOAD_FUNC(conv);
LOAD_FUNC(conv_cmplx);
LOAD_FUNC(conv_fft);
diff --git a/include/dspl.h b/include/dspl.h
index 62e4d37..13062d8 100644
--- a/include/dspl.h
+++ b/include/dspl.h
@@ -87,6 +87,27 @@ typedef double complex_t[2];
+/* Point 2D point2d_t[0] - x
+ point2d_t[1] - y
+*/
+typedef double point2d_t[2];
+
+typedef struct
+{
+ point2d_t* points; /* line points array */
+ int npoints; /* number of points */
+}line2d_t;
+
+typedef struct
+{
+ line2d_t* lines; /* lines array */
+ int nlines; /* number of lines */
+ double level; /* contour level */
+}contour2d_t;
+
+
+
+
#ifdef DOXYGEN_ENGLISH
/*! ****************************************************************************
\ingroup DFT_GROUP
@@ -836,6 +857,14 @@ DECLARE_FUNC(int, concat, void*
COMMA size_t
COMMA void*);
/*----------------------------------------------------------------------------*/
+DECLARE_FUNC(int, contour2d, double* z
+ COMMA double* x
+ COMMA double* y
+ COMMA int n
+ COMMA int m
+ COMMA double lev
+ COMMA contour2d_t* c);
+/*----------------------------------------------------------------------------*/
DECLARE_FUNC(int, conv, double*
COMMA int
COMMA double*