From 7773d2867d3287581cf55284f572dd79c71c9dd7 Mon Sep 17 00:00:00 2001 From: Dsplib Date: Mon, 8 Jan 2024 22:57:31 +0300 Subject: [PATCH] new file: dspl/src/gnuplot/contour2d.c modified: examples/src/contour_test.c modified: include/dspl.c modified: include/dspl.h --- dspl/src/gnuplot/contour2d.c | 558 +++++++++++++++++++++++++++++++++++ examples/src/contour_test.c | 374 +---------------------- include/dspl.c | 2 + include/dspl.h | 29 ++ 4 files changed, 594 insertions(+), 369 deletions(-) create mode 100644 dspl/src/gnuplot/contour2d.c 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*