Skip to content
Snippets Groups Projects
Commit ccdec088 authored by slevy's avatar slevy
Browse files

warp -wx a:val, b:val, ... -wy ... -wz ... -add dx,dy,dz 

Works but no -ref yet.
parent 2c45228b
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 109 additions and 79 deletions
src/warp.c
+ 109
− 79
View file @ ccdec088
......@@ -19,6 +19,7 @@
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include "shmem.h"
#include "specks.h"
......@@ -39,12 +40,18 @@ struct speckcache {
typedef enum { DIFFROT, EXTRAPOLATE, SUBSTITUTE } warpstyle_t;
struct fvvec {
float *v;
int len, room;
struct svvec {
struct cf {
int c;
float f;
} *v;
int len, room, maxc;
};
struct warpstuff {
struct stuff *st;
warpstyle_t style;
/* for style == DIFFROT (differential rotation in XZ plane) */
......@@ -72,9 +79,9 @@ struct warpstuff {
/* style == SUBSTITUTE ("-subst field0") */
/* Substitute some sequence of 3 attributes -> XYZ position */
int subst0;
struct fvvec sref, sadd;
struct fvvec sproj[3];
struct svvec sref;
struct svvec sproj[3];
struct svvec sadd;
int valid;
int locked;
......@@ -161,70 +168,96 @@ int parsematrix( double tfm[16], char *str, const char *whose )
return 1;
}
static void vvneed( struct fvvec *fv, int need )
static void svneed( struct svvec *sv, int need )
{
if(need < fv->room)
if(need < sv->room)
return;
fv->room = 3 + 2 * (need > fv->room ? need : fv->room);
fv->v = (float *) (fv->v ? malloc( fv->room * sizeof(float) )
: realloc( fv->v, fv->room * sizeof(float) ));
memset( &fv->v[fv->len], 0, (fv->room - fv->len) * sizeof(float) );
sv->room = 3 + 2 * (need > sv->room ? need : sv->room);
sv->v = (struct cf *) (sv->v ? malloc( sv->room * sizeof(*sv->v) )
: realloc( sv->v, sv->room * sizeof(*sv->v) ));
memset( &sv->v[sv->len], 0, (sv->room - sv->len) * sizeof(*sv->v) );
}
float *fvvref( struct fvvec *fv, int index )
struct cf *svvref( struct svvec *sv, int index )
{
if(index < 0)
return NULL;
if(index >= fv->room)
vvneed( fv, index+1 );
return &fv->v[index];
if(index >= sv->room)
svneed( sv, index+1 );
if(sv->len <= index)
sv->len = index+1;
return &sv->v[index];
}
static int substvec( struct warpstuff *ws, struct fvvec *fv, int arg, int argc, char **argv )
static int substvec( struct warpstuff *ws, struct svvec *sv, int arg, int argc, char **argv )
{
int i, a;
int ix, ok, a;
char *cp, *ep;
int err = 0;
int startcoef = -1;
int i = 0;
int coef = 0;
ws->style = SUBSTITUTE;
i = 0;
sv->maxc = -1;
a = arg;
/*
* expect: [integer ":"] float float ... {EOF | -<non-number>}
* expect: [integer-or-fieldname ":"] float float ... {EOF | -<non-number>}
*/
cp = argv[a];
do {
double v = strtod(cp, &ep);
if(ep == cp) {
double v;
char *colon;
struct cf *cfp;
while(*cp == ',' || isspace(*cp))
cp++;
if(*cp == '\0') {
if(++a >= argc)
break;
if(argv[a][0] == '-' && !(argv[a][1] == '.' || isdigit(argv[a][1])))
break;
err = 1;
break;
cp = argv[a];
continue;
}
if(*ep == ':') {
if(i != 0) {
msg("warp %s ... %s -- expected startindex: coef0 coef1 ... or just list of coefs. What's this?",
argv[arg-1], argv[a]);
err = 1;
colon = strpbrk(cp, ":, \t");
if(colon && *colon == ':') {
*colon = '\0';
if(colon == cp+1 && *cp>='x' && *cp<='z') {
coef = &(((struct speck *)0)->p.x[*cp - 'x']) - &(((struct speck *)0)->val[0]);
} else {
err = !specks_set_byvariable( ws->st, cp, &coef );
}
*colon = ':';
if(err) {
msg("warp %s: expected {<integer>|<fieldname>|x|y|z}: coef, coef, ..., not %s",
argv[arg-1], cp);
break;
}
startcoef = (int) v;
cp = colon+1;
continue;
}
} else {
for(cp = ep; *cp == ',' || isspace(*cp); cp++)
;
if(*cp == '\0') {
if(++a >= argc)
break;
cp = argv[a];
}
*fvvref( fv, i ) = v;
i++;
v = strtod(cp, &ep);
if(ep == cp) {
err = 1;
break;
}
cp = ep;
cfp = svvref( sv, i );
cfp->f = v;
cfp->c = coef;
if(sv->maxc < coef)
sv->maxc = coef;
coef++;
i++;
} while(!err);
return err;
return err ? argc+1 : a;
}
struct warpstuff *warp_setup( struct stuff *st, struct warpstuff *oldws, int argc, char **argv )
......@@ -241,24 +274,22 @@ struct warpstuff *warp_setup( struct stuff *st, struct warpstuff *oldws, int arg
tws = *oldws;
else
memset(&tws, 0, sizeof(tws));
tws.st = st;
for(i = 1; i < argc-1; i += 2) {
char *optarg = argv[i+1];
if(!strncmp(argv[i], "-substitute", 3)) {
char subst0name[32];
tws.style = SUBSTITUTE;
sscanf(optarg, "%31[^,]", subst0name);
if(!specks_set_byvariable( st, subst0name, &tws.subst0 )) {
msg("warp -substitute: expected fieldnumber or name, not \"%s\"", subst0name);
}
} else if(!strncmp(argv[i], "-w", 2) && argv[i][2] >= 'x' && argv[i][2] <= 'z') {
i = substvec( &tws, &tws.sproj[ argv[i][2] - 'x' ], i, argc, argv );
if(!strncmp(argv[i], "-w", 2) && argv[i][2] >= 'x' && argv[i][2] <= 'z') {
i = substvec( &tws, &tws.sproj[ argv[i][2] - 'x' ], i+1, argc, argv ) - 2;
} else if(!strcmp(argv[i], "-ref")) {
i = substvec( &tws, &tws.sref, i, argc, argv );
i = substvec( &tws, &tws.sref, i+1, argc, argv ) - 2;
} else if(!strcmp(argv[i], "-add")) {
i = substvec( &tws, &tws.sadd, i, argc, argv );
i = substvec( &tws, &tws.sadd, i+1, argc, argv ) - 2;
if(tws.sadd.len != 3 || tws.sadd.v[0].c != 0 || tws.sadd.v[1].c != 1 || tws.sadd.v[2].c != 2) {
msg("warp -add: expected dx,dy,dz");
tws.sadd.len = 0;
}
} else if(!strncmp(argv[i], "-extrapolate", 3)) {
char coef0name[32];
......@@ -467,35 +498,34 @@ void warpspecks( struct warpstuff *ws,
case SUBSTITUTE: {
int k;
int maxval = SPECKMAXVAL(osl);
int maxm = maxval - ws->subst0;
int maxref = ws->sref.len < maxm ? ws->sref.len : maxm;
for(i = 0; i < n; i++) {
osp = NextSpeck(osp, osl, 1);
sp = NextSpeck(sp, sl, 1);
for(k = 0; k < 3; k++) {
struct fvvec *sproj = &ws->sproj[k];
float *oval = &osp->val[ws->subst0];
float v = ws->sadd.len>k ? ws->sadd.v[k] : 0;
if(sproj->len == 0) {
if(k < maxm)
v += oval[k];
} else {
int m, maxcoef;
maxcoef = sproj->len < maxm ? sproj->len : maxm;
if(maxref == 0) {
for(m = 0; m < maxcoef; m++)
v += sproj->v[m] * oval[m];
} else if(maxref >= maxcoef) {
for(m = 0; m < maxcoef; m++)
v += sproj->v[m] * (oval[m] - ws->sref.v[m]);
int maxcoef = ws->sproj[0].maxc > ws->sproj[1].maxc ? ws->sproj[0].maxc : ws->sproj[1].maxc;
if(maxcoef < ws->sproj[2].maxc)
maxcoef = ws->sproj[2].maxc;
if(maxcoef >= maxval) {
/* Can't transform, or identity tfm */
for(i = 0; i < n; i++) {
sp->p = osp->p;
osp = NextSpeck(osp, osl, 1);
sp = NextSpeck(sp, sl, 1);
}
} else {
for(i = 0; i < n; i++) {
for(k = 0; k < 3; k++) {
struct svvec *sproj = &ws->sproj[k];
float *oval = &osp->val[0];
float v = ws->sadd.len ? ws->sadd.v[k].f : 0; /* we guarantee sadd.len == 0 or 3 */
if(sproj->len == 0) {
v += osp->p.x[k];
} else {
for(m = 0; m < maxcoef; m++)
v += sproj->v[m] *
(m<maxref ? (oval[m] - ws->sref.v[m]) : oval[m]);
int m;
for(m = 0; m < sproj->len; m++)
v += osp->val[ sproj->v[m].c ] * sproj->v[m].f;
}
sp->p.x[k] = v;
}
sp->p.x[k] = v;
osp = NextSpeck(osp, osl, 1);
sp = NextSpeck(sp, sl, 1);
}
}
}
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

Contact help@ncsa.illinois.edu with questions regarding this site. | All rights reserved. ©2024 Board of Trustees of the University of Illinois. | Web Privacy Notice