Изменения, которые будут включены в коммит:

изменено:      dspl/lapack/la_xisnan.mod
	изменено:      examples/src/contour_test.c
master
Dsplib 2023-12-20 23:07:39 +03:00
rodzic 6726287f55
commit 6c8cbfa557
2 zmienionych plików z 393 dodań i 339 usunięć

Plik binarny nie jest wyświetlany.

Wyświetl plik

@ -1,339 +1,393 @@
#include <stdio.h> #include <stdio.h>
#include <stdlib.h> #include <stdlib.h>
#include <string.h> #include <string.h>
#include "dspl.h" #include "dspl.h"
#define N 3 #define N 3
#define M 4 #define M 4
#define BSIZE 256 #define BSIZE 256
typedef double point2d_t[2]; typedef double point2d_t[2];
typedef struct typedef struct
{ {
point2d_t p[2]; point2d_t p[2];
int flag; int flag;
} linseg_t; } linseg_t;
int add_linseg(linseg_t** ls, int* lsnum, int* lscnt, int add_linseg(linseg_t** ls, int* lsnum, int* lscnt,
point2d_t* p0, point2d_t* p1) point2d_t* p0, point2d_t* p1)
{ {
int n, c; int n, c;
n = *lsnum; n = *lsnum;
c = *lscnt; c = *lscnt;
// проверяем выделение памяти. // проверяем выделение памяти.
if((n == 0) && ((*ls)==NULL)) if((n == 0) && ((*ls)==NULL))
{ {
n = BSIZE; n = BSIZE;
(*ls) = (linseg_t*)malloc(n * sizeof(linseg_t)); (*ls) = (linseg_t*)malloc(n * sizeof(linseg_t));
} }
else else
{ {
// при необходимости увеличиваем // при необходимости увеличиваем
if(c >= n) if(c >= n)
{ {
n += BSIZE; n += BSIZE;
(*ls) = (linseg_t*)realloc((*ls), n * sizeof(linseg_t)); (*ls) = (linseg_t*)realloc((*ls), n * sizeof(linseg_t));
} }
} }
(*ls)[c].p[0][0] = p0[0][0]; (*ls)[c].p[0][0] = p0[0][0];
(*ls)[c].p[0][1] = p0[0][1]; (*ls)[c].p[0][1] = p0[0][1];
(*ls)[c].p[1][0] = p1[0][0]; (*ls)[c].p[1][0] = p1[0][0];
(*ls)[c].p[1][1] = p1[0][1]; (*ls)[c].p[1][1] = p1[0][1];
(*ls)[c].flag = 1; (*ls)[c].flag = 1;
c++; c++;
(*lsnum) = n; (*lsnum) = n;
(*lscnt) = c; (*lscnt) = c;
return RES_OK; return RES_OK;
} }
int linseg_create(double* z, double* x, double* y, int linseg_create(double* z, double* x, double* y,
int n, int m, double lev, int n, int m, double lev,
linseg_t** ls, int* sz) linseg_t** ls, int* sz)
{ {
int lsnum, lscnt, t, in, im, i; int lsnum, lscnt, t, in, im, i;
point2d_t p0 = {0}; point2d_t p0 = {0};
point2d_t p1 = {0}; point2d_t p1 = {0};
double dx; double dx;
double dy; double dy;
if((ls== NULL)||(z==NULL)) if((ls== NULL)||(z==NULL))
return ERROR_PTR; return ERROR_PTR;
lsnum = 0; lsnum = 0;
lscnt = 0; lscnt = 0;
for(in = 0; in < n-1; in++) for(in = 0; in < n-1; in++)
{ {
for(im = 0; im < m-1; im++) for(im = 0; im < m-1; im++)
{ {
i = in + im * n; i = in + im * n;
t = 0; t = 0;
t += z[i] > lev ? 8 : 0; t += z[i] > lev ? 8 : 0;
t += z[i+n] > lev ? 4 : 0; t += z[i+n] > lev ? 4 : 0;
t += z[i+n+1] > lev ? 2 : 0; t += z[i+n+1] > lev ? 2 : 0;
t += z[i+1] > lev ? 1 : 0; t += z[i+1] > lev ? 1 : 0;
printf("%d, %d, %d\n", in, im, t); printf("%d, %d, %d\n", in, im, t);
switch(t) switch(t)
{ {
case 0: case 0:
case 15: case 15:
break; break;
case 1: case 1:
case 14: case 14:
// z[i,j] * (1-dx) + z[i+1, j] * dx = lev // z[i,j] * (1-dx) + z[i+1, j] * dx = lev
// z[i,j] + dx *(z[i+1, j] - z[i,j]) = lev // z[i,j] + dx *(z[i+1, j] - z[i,j]) = lev
dx = (lev - z[i]) / (z[i+1] - z[i]); dx = (lev - z[i]) / (z[i+1] - z[i]);
p0[0] = x[in] + dx; p0[0] = x[in] + dx;
p0[1] = y[im]; p0[1] = y[im];
dy = (lev - z[i+1]) / (z[i+n+1] - z[i+1]); dy = (lev - z[i+1]) / (z[i+n+1] - z[i+1]);
p1[0] = x[in+1]; p1[0] = x[in+1];
p1[1] = y[im] + dy; p1[1] = y[im] + dy;
add_linseg(ls, &lsnum, &lscnt, &p0, &p1); add_linseg(ls, &lsnum, &lscnt, &p0, &p1);
break; break;
case 2: case 2:
case 13: case 13:
dx = (lev - z[i+n]) / (z[i+n+1] - z[i+n]); dx = (lev - z[i+n]) / (z[i+n+1] - z[i+n]);
p0[0] = x[in] + dx; p0[0] = x[in] + dx;
p0[1] = y[im + 1]; p0[1] = y[im + 1];
dy = (lev - z[i+1]) / (z[i+n+1] - z[i+1]); dy = (lev - z[i+1]) / (z[i+n+1] - z[i+1]);
p1[0] = x[in + 1]; p1[0] = x[in + 1];
p1[1] = y[im] + dy; p1[1] = y[im] + dy;
add_linseg(ls, &lsnum, &lscnt, &p0, &p1); add_linseg(ls, &lsnum, &lscnt, &p0, &p1);
break; break;
case 3: case 3:
case 12: case 12:
dx = (lev - z[i]) / (z[i+1] - z[i]); dx = (lev - z[i]) / (z[i+1] - z[i]);
p0[0] = x[in] + dx; p0[0] = x[in] + dx;
p0[1] = y[im]; p0[1] = y[im];
dx = (lev - z[i+n]) / (z[i+n+1] - z[i+n]); dx = (lev - z[i+n]) / (z[i+n+1] - z[i+n]);
p1[0] = x[in] + dx; p1[0] = x[in] + dx;
p1[1] = y[im+1]; p1[1] = y[im+1];
add_linseg(ls, &lsnum, &lscnt, &p0, &p1); add_linseg(ls, &lsnum, &lscnt, &p0, &p1);
break; break;
case 4: case 4:
case 11: case 11:
dy = (lev - z[i]) / (z[i+n] - z[i]); dy = (lev - z[i]) / (z[i+n] - z[i]);
p0[0] = x[in]; p0[0] = x[in];
p0[1] = y[im] + dy; p0[1] = y[im] + dy;
dx = (lev - z[i+n]) / (z[i+n+1] - z[i+n]); dx = (lev - z[i+n]) / (z[i+n+1] - z[i+n]);
p1[0] = x[in] + dx; p1[0] = x[in] + dx;
p1[1] = y[im + 1]; p1[1] = y[im + 1];
add_linseg(ls, &lsnum, &lscnt, &p0, &p1); add_linseg(ls, &lsnum, &lscnt, &p0, &p1);
break; break;
case 5: case 5:
dy = (lev - z[i]) / (z[i+n] - z[i]); dy = (lev - z[i]) / (z[i+n] - z[i]);
p0[0] = x[in]; p0[0] = x[in];
p0[1] = y[im] + dy; p0[1] = y[im] + dy;
dx = (lev - z[i]) / (z[i+1] - z[i]); dx = (lev - z[i]) / (z[i+1] - z[i]);
p1[0] = x[in] + dx; p1[0] = x[in] + dx;
p1[1] = y[im]; p1[1] = y[im];
add_linseg(ls, &lsnum, &lscnt, &p0, &p1); add_linseg(ls, &lsnum, &lscnt, &p0, &p1);
dx = (lev - z[i+n]) / (z[i+n+1] - z[i+n]); dx = (lev - z[i+n]) / (z[i+n+1] - z[i+n]);
p0[0] = x[in]+dx; p0[0] = x[in]+dx;
p0[1] = y[im+1]; p0[1] = y[im+1];
dy = (lev - z[i+1]) / (z[i+n+1] - z[i+1]); dy = (lev - z[i+1]) / (z[i+n+1] - z[i+1]);
p1[0] = x[in+1]; p1[0] = x[in+1];
p1[1] = y[im] + dy; p1[1] = y[im] + dy;
add_linseg(ls, &lsnum, &lscnt, &p0, &p1); add_linseg(ls, &lsnum, &lscnt, &p0, &p1);
break; break;
case 6: case 6:
case 9: case 9:
dy = (lev - z[i]) / (z[i+n] - z[i]); dy = (lev - z[i]) / (z[i+n] - z[i]);
p0[0] = x[in]; p0[0] = x[in];
p0[1] = y[im] + dy; p0[1] = y[im] + dy;
dy = (lev - z[i+1]) / (z[i+n+1] - z[i+1]); dy = (lev - z[i+1]) / (z[i+n+1] - z[i+1]);
p1[0] = x[in+1]; p1[0] = x[in+1];
p1[1] = y[im]; p1[1] = y[im]+dy;
add_linseg(ls, &lsnum, &lscnt, &p0, &p1); add_linseg(ls, &lsnum, &lscnt, &p0, &p1);
break; break;
case 7: case 7:
case 8: case 8:
dy = (lev - z[i]) / (z[i+n] - z[i]); dy = (lev - z[i]) / (z[i+n] - z[i]);
p0[0] = x[in]; p0[0] = x[in];
p0[1] = y[im] + dy; p0[1] = y[im] + dy;
dx = (lev - z[i]) / (z[i+1] - z[i]); dx = (lev - z[i]) / (z[i+1] - z[i]);
p1[0] = x[in] + dx; p1[0] = x[in] + dx;
p1[1] = y[im]; p1[1] = y[im];
add_linseg(ls, &lsnum, &lscnt, &p0, &p1); add_linseg(ls, &lsnum, &lscnt, &p0, &p1);
break; break;
case 10: case 10:
dy = (lev - z[i]) / (z[i+n] - z[i]); dy = (lev - z[i]) / (z[i+n] - z[i]);
p0[0] = x[in]; p0[0] = x[in];
p0[1] = y[im] + dy; p0[1] = y[im] + dy;
dx = (lev - z[i+n]) / (z[i+n+1] - z[i+1]); dx = (lev - z[i+n]) / (z[i+n+1] - z[i+1]);
p1[0] = x[in] + dx; p1[0] = x[in] + dx;
p1[1] = y[im+1]; p1[1] = y[im+1];
add_linseg(ls, &lsnum, &lscnt, &p0, &p1); add_linseg(ls, &lsnum, &lscnt, &p0, &p1);
dx = (lev - z[i]) / (z[i+1] - z[i]); dx = (lev - z[i]) / (z[i+1] - z[i]);
p0[0] = x[in]+dx; p0[0] = x[in]+dx;
p0[1] = y[im]; p0[1] = y[im];
dy = (lev - z[i+1]) / (z[i+n+1] - z[i+1]); dy = (lev - z[i+1]) / (z[i+n+1] - z[i+1]);
p1[0] = x[in+1]; p1[0] = x[in+1];
p1[1] = y[im] + dy; p1[1] = y[im] + dy;
add_linseg(ls, &lsnum, &lscnt, &p0, &p1); add_linseg(ls, &lsnum, &lscnt, &p0, &p1);
break; break;
default: default:
break; break;
} }
} }
} }
*ls = (linseg_t*)realloc(*ls, lscnt * sizeof(linseg_t)); *ls = (linseg_t*)realloc(*ls, lscnt * sizeof(linseg_t));
*sz = lscnt; *sz = lscnt;
return RES_OK; return RES_OK;
} }
double dist(point2d_t* p0, point2d_t* p1) double dist(point2d_t* p0, point2d_t* p1)
{ {
double dx, dy; double dx, dy;
dx = p0[0][0] - p1[0][0]; dx = p0[0][0] - p1[0][0];
dy = p0[0][1] - p1[0][1]; dy = p0[0][1] - p1[0][1];
return sqrt(dx*dx + dy*dy); return sqrt(dx*dx + dy*dy);
} }
int line_create(linseg_t* ls, int nls, point2d_t** line, int* np) int line_create(linseg_t* ls, int nls, point2d_t** line, int* np)
{ {
int i, j, c, n; int i, j, c, n, k;
if(!line || !ls || !np) if(!line || !ls || !np)
return ERROR_PTR; return ERROR_PTR;
// printf("*line = %x, *np = %d \n", *line, *np); // printf("*line = %x, *np = %d \n", *line, *np);
i = 0; i = 0;
while(!(ls[i].flag) && i < nls) while(!(ls[i].flag) && i < nls)
i++; i++;
if(i==nls) if(i==nls)
{ {
*np = 0; *np = 0;
return RES_OK; return RES_OK;
} }
//printf("i = %d, ls[%d] = [%.1f %.1f] -- [%.1f, %.1f]\n", //printf("i = %d, ls[%d] = [%.1f %.1f] -- [%.1f, %.1f]\n",
// i, i, ls[i].p[0][0], ls[i].p[0][1], // i, i, ls[i].p[0][0], ls[i].p[0][1],
// ls[i].p[1][0], ls[i].p[1][1]); // ls[i].p[1][0], ls[i].p[1][1]);
n = BSIZE; n = BSIZE;
if((*line) == NULL) if((*line) == NULL)
(*line) = (point2d_t*)malloc(n*sizeof(point2d_t)); (*line) = (point2d_t*)malloc(n*sizeof(point2d_t));
else else
(*line) = (point2d_t*)realloc((*line), n*sizeof(point2d_t)); (*line) = (point2d_t*)realloc((*line), n*sizeof(point2d_t));
c = 0; c = 0;
(*line)[c][0] = ls[i].p[0][0]; (*line)[c][0] = ls[i].p[0][0];
(*line)[c][1] = ls[i].p[0][1]; (*line)[c][1] = ls[i].p[0][1];
c++; c++;
(*line)[c][0] = ls[i].p[1][0]; (*line)[c][0] = ls[i].p[1][0];
(*line)[c][1] = ls[i].p[1][1]; (*line)[c][1] = ls[i].p[1][1];
c++; c++;
ls[i].flag = 0; ls[i].flag = 0;
for(i = 0; i < nls; i++) for(i = 0; i < nls; i++)
{ {
for(j = 0; j < nls; j++) for(j = 0; j < nls; j++)
{ {
if(ls[j].flag) if(ls[j].flag)
{ {
//сравниваем с первой точкой отрезка ls[j] //сравниваем последнюю точку линии с первой точкой отрезка ls[j]
if(dist((*line)+c-1, ls[j].p) < 1E-8) if(dist((*line)+c-1, ls[j].p) < 1E-8)
{ {
// printf("c0 = %d, j0 = %d\n", c, j);
// проверяем выделение памяти.
// проверяем выделение памяти. // при необходимости увеличиваем
// при необходимости увеличиваем if(c>=n)
if(c>=n) {
{ n += BSIZE;
n += BSIZE; (*line) = (point2d_t*)realloc(*line, n*sizeof(point2d_t));
(*line) = (point2d_t*)realloc(*line, n*sizeof(point2d_t)); }
} //если первая точка отрезка совпадает, то добавляем
//если первая точка совпадает, то добавляем //в линию вторую точку СЗАДИ
//в линию вторую точку (*line)[c][0] = ls[j].p[1][0];
(*line)[c][0] = ls[j].p[1][0]; (*line)[c][1] = ls[j].p[1][1];
(*line)[c][1] = ls[j].p[1][1]; ls[j].flag = 0;
ls[j].flag = 0; c++;
c++; }
} }
} if(ls[j].flag)
if(ls[j].flag) {
{ // сравниваем последнюю точку линии
//сравниваем со второй точкой отрезка ls[j] // со второй точкой отрезка ls[j]
if(dist((*line)+c-1, ls[j].p+1) < 1E-8) if(dist((*line)+c-1, ls[j].p+1) < 1E-8)
{ {
// printf("c1 = %d, j1 = %d\n", c, j); // проверяем выделение памяти.
// при необходимости увеличиваем
// проверяем выделение памяти. if(c>=n)
// при необходимости увеличиваем {
if(c>=n) n += BSIZE;
{ (*line) = (point2d_t*)realloc(*line, n*sizeof(point2d_t));
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];
(*line)[c][0] = ls[j].p[0][0]; ls[j].flag = 0;
(*line)[c][1] = ls[j].p[0][1]; c++;
ls[j].flag = 0; }
c++; }
}
} if(ls[j].flag)
} {
} //сравниваем последнюю точку линии с первой точкой отрезка ls[j]
(*line) = (point2d_t*)realloc((*line), c*sizeof(point2d_t)); if(dist((*line), ls[j].p) < 1E-8)
*np = c; {
return RES_OK; // проверяем выделение памяти.
} // при необходимости увеличиваем
if(c>=n)
{
n += BSIZE;
(*line) = (point2d_t*)realloc(*line, n*sizeof(point2d_t));
}
int main(int argc, char* argv[]) //если первая точка отрезка совпадает, то добавляем
{ //в линию вторую точку СПЕРЕДИ
/* Matrix for(k = c; k>0; k--)
z = [ 0 0 1 0; {
0 1 1 0; (*line)[k][0] = (*line)[k-1][0];
0 0 0 0]; (*line)[k][1] = (*line)[k-1][1];
in array a by columns }
*/ (*line)[0][0] = ls[j].p[1][0];
double z[N*M] = { 0, 0, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0}; (*line)[0][1] = ls[j].p[1][1];
double x[N] = {0.0, 1.0, 2.0}; ls[j].flag = 0;
double y[M] = {0.0, 1.0, 2.0, 3.0}; c++;
linseg_t *ls = NULL; }
int nls, i, j; }
point2d_t** line = NULL; if(ls[j].flag)
int *np = NULL; {
// сравниваем последнюю точку линии
line = (point2d_t**)malloc(10*sizeof(point2d_t*)); // со второй точкой отрезка ls[j]
memset(line,0,10*sizeof(point2d_t*)); if(dist((*line), ls[j].p+1) < 1E-8)
{
np = (int*)malloc(10*sizeof(int)); // проверяем выделение памяти.
memset(np, 0, 10*sizeof(int)); // при необходимости увеличиваем
if(c>=n)
linseg_create(z, x, y, N, M, 0.5, &ls, &nls); {
for(i =0; i< nls; i++) n += BSIZE;
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]); (*line) = (point2d_t*)realloc(*line, n*sizeof(point2d_t));
}
//если вторая точка совпадает, то добавляем
i = -1; //в линию первую точку СПЕРЕДИ
do
{ for(k = c; k>0; k--)
i++; {
line_create(ls, nls, line+i, np); (*line)[k][0] = (*line)[k-1][0];
for(j =0; j< np[i]; j++) (*line)[k][1] = (*line)[k-1][1];
printf("%[%.1f %.1f] -- ", line[i][j][0], line[i][j][1]); }
printf("\n");
}while(np[i]); (*line)[0][0] = ls[j].p[0][0];
(*line)[0][1] = ls[j].p[0][1];
ls[j].flag = 0;
return 0; c++;
} }
}
}
}
(*line) = (point2d_t*)realloc((*line), c*sizeof(point2d_t));
*np = c;
return RES_OK;
}
int main(int argc, char* argv[])
{
/* Matrix
z = [ 0 0 1 1;
0 1 0 0;
0 1 0 0];
in array a by columns
*/
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));
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]);
return 0;
}