#include "moss.h" void mossProjTran(float *xOutP, float *yOutP, float m[3][3], float xIn, float yIn) { float t[3]; t[0] = m[0][0]*xIn + m[0][1]*yIn + m[0][2]; t[1] = m[1][0]*xIn + m[1][1]*yIn + m[1][2]; t[2] = m[2][0]*xIn + m[2][1]*yIn + m[2][2]; *xOutP = t[0]/t[2]; *yOutP = t[1]/t[2]; } void mossImageSize(int *sizeXP, int *sizeYP, Nrrd *img) { int base; if (2 == img->dim) { base = 0; } else { base = 1; } *sizeXP = img->size[0+base]; *sizeYP = img->size[1+base]; } void mossImageBounds(float *minXP, float *maxXP, float *minYP, float *maxYP, Nrrd *img) { int base; if (2 == img->dim) { base = 0; } else { base = 1; } *minXP = img->axisMin[0+base]; *maxXP = img->axisMax[0+base]; *minYP = img->axisMin[1+base]; *maxYP = img->axisMax[1+base]; } int mossInterp(float *wght, float *val, Nrrd *img, float xs, float ys) { int sx, sy, xi, yi; unsigned char *data; float xf, yf, w00, w01, w10, w11, wx, wy; data = (unsigned char *)img->data; if (2 == img->dim) { sx = img->size[0]; sy = img->size[1]; } else { sx = img->size[1]; sy = img->size[2]; } if (!(AIR_INSIDE(0, xs, sx-1) && AIR_INSIDE(0, ys, sy-1))) { wght[0] = 0; val[0] = val[1] = val[2] = 0; return 0; } xi = xs; yi = ys; xi -= (xi == sx-1); xf = xs - xi; yi -= (yi == sy-1); yf = ys - yi; w00 = (1-xf)*(1-yf); w01 = xf*(1-yf); w10 = (1-xf)*yf; w11 = xf*yf; if (2 == img->dim) { val[0] = (w00*data[xi+0 + sx*(yi+0)] + w01*data[xi+1 + sx*(yi+0)] + w10*data[xi+0 + sx*(yi+1)] + w11*data[xi+1 + sx*(yi+1)]); } else { val[0] = (w00*data[0 + 3*(xi+0 + sx*(yi+0))] + w01*data[0 + 3*(xi+1 + sx*(yi+0))] + w10*data[0 + 3*(xi+0 + sx*(yi+1))] + w11*data[0 + 3*(xi+1 + sx*(yi+1))]); val[1] = (w00*data[1 + 3*(xi+0 + sx*(yi+0))] + w01*data[1 + 3*(xi+1 + sx*(yi+0))] + w10*data[1 + 3*(xi+0 + sx*(yi+1))] + w11*data[1 + 3*(xi+1 + sx*(yi+1))]); val[2] = (w00*data[2 + 3*(xi+0 + sx*(yi+0))] + w01*data[2 + 3*(xi+1 + sx*(yi+0))] + w10*data[2 + 3*(xi+0 + sx*(yi+1))] + w11*data[2 + 3*(xi+1 + sx*(yi+1))]); } wx = AIR_AFFINE(0, xs, sx-1, 0, 2); wx = AIR_CLAMP(0, wx, 2); if (wx > 1) wx = 2 - wx; wy = AIR_AFFINE(0, ys, sy-1, 0, 2); wy = AIR_CLAMP(0, wy, 2); if (wy > 1) wy = 2 - wy; *wght = AIR_MIN(wx, wy); /* wght goes from 0.0 to 1.0 linearly, we want to send it through a bit of a sinusoid */ *wght = sin(M_PI*(*wght)/2); *wght *= *wght; if (*wght) return 1; else return 0; } void mossSetVal(Nrrd *img, int x, int y, float v0, float v1, float v2) { int sx; unsigned char *data; data = (unsigned char *)img->data; if (2 == img->dim) { sx = img->size[0]; data[x + sx*y] = v0; } else { sx = img->size[1]; data[0 + 3*(x + sx*y)] = v0; data[1 + 3*(x + sx*y)] = v1; data[2 + 3*(x + sx*y)] = v2; } } int mossLearnBounds(float *minXP, float *maxXP, float *minYP, float *maxYP, mossMatrix *matx, Nrrd **imgs, int numImgs, int which) { char me[]="mossLearnBounds", err[128]; int i, sizeX, sizeY; float x, y, (*m)[3][3]; mossCorresp *pair; if (!(minXP && maxXP && minYP && maxYP && matx && imgs)) { sprintf(err, BIFF_NULL, me); biffAdd(MOSS, err); return 1; } *minXP = *minYP = airPosInff(); *maxXP = *maxYP = airNegInff(); for (i=0; i<=numImgs-1; i++) { pair = matx->pair[i][which]; if (i < which) m = &(pair->m01); else m = &(pair->m10); mossImageSize(&sizeX, &sizeY, imgs[i]); printf("%s: image %d has size %d %d\n", me, i, sizeX, sizeY); mossProjTran(&x, &y, *m, 0, 0); *minXP = AIR_MIN(*minXP, x); *minYP = AIR_MIN(*minYP, y); *maxXP = AIR_MAX(*maxXP, x); *maxYP = AIR_MAX(*maxYP, y); mossProjTran(&x, &y, *m, sizeX-1, 0); *minXP = AIR_MIN(*minXP, x); *minYP = AIR_MIN(*minYP, y); *maxXP = AIR_MAX(*maxXP, x); *maxYP = AIR_MAX(*maxYP, y); mossProjTran(&x, &y, *m, sizeX-1, sizeY-1); *minXP = AIR_MIN(*minXP, x); *minYP = AIR_MIN(*minYP, y); *maxXP = AIR_MAX(*maxXP, x); *maxYP = AIR_MAX(*maxYP, y); mossProjTran(&x, &y, *m, 0, sizeY-1); *minXP = AIR_MIN(*minXP, x); *minYP = AIR_MIN(*minYP, y); *maxXP = AIR_MAX(*maxXP, x); *maxYP = AIR_MAX(*maxYP, y); } return 0; } int mossFillImage(Nrrd *out, mossMatrix *matx, Nrrd **imgs, int numImgs, int which) { char me[]="mossFillImage", err[128]; int idx, incr, i, x, y, sx, sy, inside; float minX, minY, maxX, maxY, xt, yt, xs, ys, *wght, (*m)[3][3], *val, totalW, totalV[3]; if (!(out && matx && imgs)) { sprintf(err, BIFF_NULL, me); biffAdd(MOSS, err); return 1; } mossImageSize(&sx, &sy, out); mossImageBounds(&minX, &maxX, &minY, &maxY, out); printf("%s: bounds: [%g,%g] [%g,%g]\n", me, minX, maxX, minY, maxY); wght = (float *)malloc(numImgs*sizeof(float)); val = (float *)malloc(3*numImgs*sizeof(float)); if (!(wght && val)) { sprintf(err, "%s: couldn't alloc wght and val arrays", me); biffAdd(MOSS, err); return 1; } for (y=0; y<=sy-1; y++) { yt = AIR_AFFINE(0, y, sy-1, minY, maxY); for (x=0; x<=sx-1; x++) { xt = AIR_AFFINE(0, x, sx-1, minX, maxX); inside = 0; for (i=0; i<=numImgs-1; i++) { if (i < which) m = &(matx->pair[i][which]->m10); else m = &(matx->pair[i][which]->m01); mossProjTran(&xs, &ys, *m, xt, yt); incr = mossInterp(wght + i, val + 3*i, imgs[i], xs, ys); if (incr) { inside += incr; idx = i; } } if (x == 168 && y == 240) { printf("wght: %g %g\n", wght[0], wght[1]); printf("inside = %d\n", inside); } if (!inside) { mossSetVal(out, x, y, 0, 0, 0); } else if (1 == inside) { mossSetVal(out, x, y, val[0 + 3*idx], val[1 + 3*idx], val[2 + 3*idx]); } else { /* we have multiple images to blend */ totalW = totalV[0] = totalV[1] = totalV[2] = 0; for (i=0; i<=numImgs-1; i++) { totalW += wght[i]; totalV[0] += wght[i]*val[0 + 3*i]; totalV[1] += wght[i]*val[1 + 3*i]; totalV[2] += wght[i]*val[2 + 3*i]; if (x == 168 && y == 240) { printf("totalW = %g\n", totalW); printf("totalV = %g %g %g\n", totalV[0], totalV[1], totalV[2]); } } mossSetVal(out, x, y, totalV[0]/totalW, totalV[1]/totalW, totalV[2]/totalW); } } } free(wght); free(val); return 0; }