1 |
|
|
/* |
2 |
|
|
Teem: Tools to process and visualize scientific data and images . |
3 |
|
|
Copyright (C) 2013, 2012, 2011, 2010, 2009 University of Chicago |
4 |
|
|
Copyright (C) 2008, 2007, 2006, 2005 Gordon Kindlmann |
5 |
|
|
Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998 University of Utah |
6 |
|
|
|
7 |
|
|
This library is free software; you can redistribute it and/or |
8 |
|
|
modify it under the terms of the GNU Lesser General Public License |
9 |
|
|
(LGPL) as published by the Free Software Foundation; either |
10 |
|
|
version 2.1 of the License, or (at your option) any later version. |
11 |
|
|
The terms of redistributing and/or modifying this software also |
12 |
|
|
include exceptions to the LGPL that facilitate static linking. |
13 |
|
|
|
14 |
|
|
This library is distributed in the hope that it will be useful, |
15 |
|
|
but WITHOUT ANY WARRANTY; without even the implied warranty of |
16 |
|
|
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
17 |
|
|
Lesser General Public License for more details. |
18 |
|
|
|
19 |
|
|
You should have received a copy of the GNU Lesser General Public License |
20 |
|
|
along with this library; if not, write to Free Software Foundation, Inc., |
21 |
|
|
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA |
22 |
|
|
*/ |
23 |
|
|
|
24 |
|
|
#include "moss.h" |
25 |
|
|
#include "privateMoss.h" |
26 |
|
|
|
27 |
|
|
int |
28 |
|
|
mossSamplerImageSet (mossSampler *smplr, Nrrd *image, float *bg) { |
29 |
|
|
static const char me[]="mossSamplerImageSet"; |
30 |
|
|
int ci, ncol; |
31 |
|
|
|
32 |
|
|
if (!(smplr && image)) { |
33 |
|
|
biffAddf(MOSS, "%s: got NULL pointer", me); |
34 |
|
|
return 1; |
35 |
|
|
} |
36 |
|
|
if (mossImageCheck(image)) { |
37 |
|
|
biffAddf(MOSS, "%s: ", me); |
38 |
|
|
return 1; |
39 |
|
|
} |
40 |
|
|
smplr->image = image; |
41 |
|
|
smplr->flag[mossFlagImage] = AIR_TRUE; |
42 |
|
|
ncol = MOSS_NCOL(image); |
43 |
|
|
smplr->bg = (float *)airFree(smplr->bg); |
44 |
|
|
if (bg) { |
45 |
|
|
smplr->bg = (float*)calloc(ncol, sizeof(float)); |
46 |
|
|
for (ci=0; ci<ncol; ci++) { |
47 |
|
|
smplr->bg[ci] = bg[ci]; |
48 |
|
|
} |
49 |
|
|
} |
50 |
|
|
return 0; |
51 |
|
|
} |
52 |
|
|
|
53 |
|
|
int |
54 |
|
|
mossSamplerKernelSet (mossSampler *smplr, |
55 |
|
|
const NrrdKernel *kernel, double *kparm) { |
56 |
|
|
static const char me[]="mossSamplerKernelSet"; |
57 |
|
|
unsigned int ki; |
58 |
|
|
|
59 |
|
|
if (!(smplr && kernel && kparm)) { |
60 |
|
|
biffAddf(MOSS, "%s: got NULL pointer", me); |
61 |
|
|
return 1; |
62 |
|
|
} |
63 |
|
|
smplr->kernel = kernel; |
64 |
|
|
for (ki=0; ki<kernel->numParm; ki++) { |
65 |
|
|
smplr->kparm[ki] = kparm[ki]; |
66 |
|
|
} |
67 |
|
|
smplr->flag[mossFlagKernel] = AIR_TRUE; |
68 |
|
|
return 0; |
69 |
|
|
} |
70 |
|
|
|
71 |
|
|
int |
72 |
|
|
mossSamplerUpdate (mossSampler *smplr) { |
73 |
|
|
static const char me[]="mossSamplerUpdate"; |
74 |
|
|
int ncol=0, fdiam=0; |
75 |
|
|
|
76 |
|
|
if (!(smplr)) { |
77 |
|
|
biffAddf(MOSS, "%s: got NULL pointer", me); |
78 |
|
|
return 1; |
79 |
|
|
} |
80 |
|
|
|
81 |
|
|
if (smplr->flag[mossFlagImage]) { |
82 |
|
|
ncol = MOSS_NCOL(smplr->image); |
83 |
|
|
if (ncol != smplr->ncol) { |
84 |
|
|
mossSamplerEmpty(smplr); |
85 |
|
|
smplr->ncol = ncol; |
86 |
|
|
} |
87 |
|
|
} |
88 |
|
|
if (smplr->flag[mossFlagKernel]) { |
89 |
|
|
fdiam = 2*AIR_ROUNDUP(smplr->kernel->support(smplr->kparm)); |
90 |
|
|
if (fdiam != smplr->fdiam) { |
91 |
|
|
mossSamplerEmpty(smplr); |
92 |
|
|
smplr->fdiam = fdiam; |
93 |
|
|
} |
94 |
|
|
} |
95 |
|
|
if (!(smplr->ivc)) { |
96 |
|
|
if (mossSamplerFill(smplr, fdiam, ncol)) { |
97 |
|
|
biffAddf(MOSS, "%s: ", me); |
98 |
|
|
return 1; |
99 |
|
|
} |
100 |
|
|
} |
101 |
|
|
if (nrrdBoundaryPad == smplr->boundary && !smplr->bg) { |
102 |
|
|
biffAddf(MOSS, "%s: want %s boundary behavior, but bg vector is NULL", |
103 |
|
|
me, airEnumStr(nrrdBoundary, nrrdBoundaryPad)); |
104 |
|
|
return 1; |
105 |
|
|
} |
106 |
|
|
|
107 |
|
|
return 0; |
108 |
|
|
} |
109 |
|
|
|
110 |
|
|
int |
111 |
|
|
mossSamplerSample (float *val, mossSampler *smplr, double xPos, double yPos) { |
112 |
|
|
static const char me[]="mossSamplerSample"; |
113 |
|
|
int i, xi, yi, ci, sx, sy, fdiam, frad, ncol; |
114 |
|
|
double xf, yf, tmp; |
115 |
|
|
float (*lup)(const void *v, size_t I); |
116 |
|
|
|
117 |
|
|
if (!(val && smplr)) { |
118 |
|
|
biffAddf(MOSS, "%s: got NULL pointer", me); |
119 |
|
|
return 1; |
120 |
|
|
} |
121 |
|
|
if (!(smplr->ivc)) { |
122 |
|
|
biffAddf(MOSS, "%s: given sampler not ready (no caches)", me); |
123 |
|
|
return 1; |
124 |
|
|
} |
125 |
|
|
|
126 |
|
|
/* set {x,y}Idx, set {x,y}Fslw to sample locations */ |
127 |
|
|
if (mossVerbose) { |
128 |
|
|
fprintf(stderr, "%s: pos = %g %g\n", me, xPos, yPos); |
129 |
|
|
} |
130 |
|
|
sx = MOSS_SX(smplr->image); |
131 |
|
|
sy = MOSS_SY(smplr->image); |
132 |
|
|
xi = (int)floor(xPos); xf = xPos - xi; |
133 |
|
|
yi = (int)floor(yPos); yf = yPos - yi; |
134 |
|
|
fdiam = smplr->fdiam; |
135 |
|
|
frad = fdiam/2; |
136 |
|
|
for (i=0; i<fdiam; i++) { |
137 |
|
|
smplr->xIdx[i] = xi + i - frad + 1; |
138 |
|
|
smplr->yIdx[i] = yi + i - frad + 1; |
139 |
|
|
smplr->xFslw[i] = xf - i + frad - 1; |
140 |
|
|
smplr->yFslw[i] = yf - i + frad - 1; |
141 |
|
|
} |
142 |
|
|
if (mossVerbose) { |
143 |
|
|
fprintf(stderr, " --> xIdx: %d %d ; xFsl %g %g\n", |
144 |
|
|
smplr->xIdx[0], smplr->xIdx[1], |
145 |
|
|
smplr->xFslw[0], smplr->xFslw[1]); |
146 |
|
|
fprintf(stderr, " yIdx: %d %d ; yFsl %g %g\n", |
147 |
|
|
smplr->yIdx[0], smplr->yIdx[1], |
148 |
|
|
smplr->yFslw[0], smplr->yFslw[1]); |
149 |
|
|
} |
150 |
|
|
switch(smplr->boundary) { |
151 |
|
|
case nrrdBoundaryBleed: |
152 |
|
|
for (i=0; i<fdiam; i++) { |
153 |
|
|
smplr->xIdx[i] = AIR_CLAMP(0, smplr->xIdx[i], sx-1); |
154 |
|
|
smplr->yIdx[i] = AIR_CLAMP(0, smplr->yIdx[i], sy-1); |
155 |
|
|
} |
156 |
|
|
break; |
157 |
|
|
case nrrdBoundaryWrap: |
158 |
|
|
for (i=0; i<fdiam; i++) { |
159 |
|
|
smplr->xIdx[i] = AIR_MOD(smplr->xIdx[i], sx); |
160 |
|
|
smplr->yIdx[i] = AIR_MOD(smplr->yIdx[i], sy); |
161 |
|
|
} |
162 |
|
|
break; |
163 |
|
|
case nrrdBoundaryPad: |
164 |
|
|
/* this is handled later */ |
165 |
|
|
break; |
166 |
|
|
default: |
167 |
|
|
biffAddf(MOSS, "%s: sorry, %s boundary not implemented", me, |
168 |
|
|
airEnumStr(nrrdBoundary, smplr->boundary)); |
169 |
|
|
return 1; |
170 |
|
|
} |
171 |
|
|
if (mossVerbose) { |
172 |
|
|
fprintf(stderr, " --> xIdx: %d %d ; xFsl %g %g\n", |
173 |
|
|
smplr->xIdx[0], smplr->xIdx[1], |
174 |
|
|
smplr->xFslw[0], smplr->xFslw[1]); |
175 |
|
|
} |
176 |
|
|
|
177 |
|
|
/* copy values to ivc, set {x,y}Fslw to filter sample weights */ |
178 |
|
|
lup = nrrdFLookup[smplr->image->type]; |
179 |
|
|
ncol = smplr->ncol; |
180 |
|
|
if (nrrdBoundaryPad == smplr->boundary) { |
181 |
|
|
for (yi=0; yi<fdiam; yi++) { |
182 |
|
|
for (xi=0; xi<fdiam; xi++) { |
183 |
|
|
if (AIR_IN_CL(0, smplr->xIdx[xi], sx-1) |
184 |
|
|
&& AIR_IN_CL(0, smplr->yIdx[yi], sy-1)) { |
185 |
|
|
for (ci=0; ci<ncol; ci++) { |
186 |
|
|
smplr->ivc[xi + fdiam*(yi + fdiam*ci)] = |
187 |
|
|
lup(smplr->image->data, |
188 |
|
|
ci + ncol*(smplr->xIdx[xi] + sx*smplr->yIdx[yi])); |
189 |
|
|
} |
190 |
|
|
} else { |
191 |
|
|
for (ci=0; ci<ncol; ci++) { |
192 |
|
|
smplr->ivc[xi + fdiam*(yi + fdiam*ci)] = smplr->bg[ci]; |
193 |
|
|
} |
194 |
|
|
} |
195 |
|
|
} |
196 |
|
|
} |
197 |
|
|
} else { |
198 |
|
|
for (yi=0; yi<fdiam; yi++) { |
199 |
|
|
for (xi=0; xi<fdiam; xi++) { |
200 |
|
|
for (ci=0; ci<ncol; ci++) { |
201 |
|
|
smplr->ivc[xi + fdiam*(yi + fdiam*ci)] = |
202 |
|
|
lup(smplr->image->data, |
203 |
|
|
ci + ncol*(smplr->xIdx[xi] + sx*smplr->yIdx[yi])); |
204 |
|
|
} |
205 |
|
|
} |
206 |
|
|
} |
207 |
|
|
} |
208 |
|
|
smplr->kernel->evalN_d(smplr->xFslw, smplr->xFslw, fdiam, smplr->kparm); |
209 |
|
|
smplr->kernel->evalN_d(smplr->yFslw, smplr->yFslw, fdiam, smplr->kparm); |
210 |
|
|
|
211 |
|
|
/* do convolution */ |
212 |
|
|
memset(val, 0, ncol*sizeof(float)); |
213 |
|
|
for (ci=0; ci<ncol; ci++) { |
214 |
|
|
for (yi=0; yi<fdiam; yi++) { |
215 |
|
|
tmp = 0; |
216 |
|
|
for (xi=0; xi<fdiam; xi++) { |
217 |
|
|
tmp += smplr->xFslw[xi]*smplr->ivc[xi + fdiam*(yi + fdiam*ci)]; |
218 |
|
|
} |
219 |
|
|
val[ci] += AIR_CAST(float, smplr->yFslw[yi]*tmp); |
220 |
|
|
} |
221 |
|
|
} |
222 |
|
|
|
223 |
|
|
return 0; |
224 |
|
|
} |