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 <teem/air.h> |
25 |
|
|
#include <teem/biff.h> |
26 |
|
|
#include <teem/hest.h> |
27 |
|
|
#include <teem/nrrd.h> |
28 |
|
|
|
29 |
|
|
static const char *overInfo = ( |
30 |
|
|
"Composites an RGBA nrrd over " |
31 |
|
|
"a background color (or image), after doing gamma correction, " |
32 |
|
|
"then quantizes to an 8-bit image. Actually, the " |
33 |
|
|
"input nrrd can have more than 4 values per pixel, " |
34 |
|
|
"but only the first four are used. If the RGBA nrrd " |
35 |
|
|
"is floating point, the values are taken at face value; " |
36 |
|
|
"if it is fixed point, the values interpreted as having " |
37 |
|
|
"been quantized (so that 8-bit RGBA images will act as " |
38 |
|
|
"you expect). When compositing with a background image, the given " |
39 |
|
|
"background image does not have to be the same size as the input " |
40 |
|
|
"image; it will be resampled (with linear interpolation) to fit. "); |
41 |
|
|
|
42 |
|
|
double |
43 |
|
|
docontrast(double val, double cfp, double cpow) { |
44 |
|
|
double v; |
45 |
|
|
|
46 |
|
|
if (val < cfp) { |
47 |
|
|
v = AIR_AFFINE(0.0, val, cfp, 0.0, 1.0); |
48 |
|
|
v = pow(v, cpow); |
49 |
|
|
val = AIR_AFFINE(0.0, v, 1.0, 0.0, cfp); |
50 |
|
|
} else { |
51 |
|
|
v = AIR_AFFINE(cfp, val, 1.0, 1.0, 0.0); |
52 |
|
|
v = pow(v, cpow); |
53 |
|
|
val = AIR_AFFINE(1.0, v, 0.0, cfp, 1.0); |
54 |
|
|
} |
55 |
|
|
return val; |
56 |
|
|
} |
57 |
|
|
|
58 |
|
|
int |
59 |
|
|
main(int argc, const char *argv[]) { |
60 |
|
|
hestOpt *hopt=NULL; |
61 |
|
|
Nrrd *nin, *nout, /* initial input and final output */ |
62 |
|
|
*ninD, /* input converted to double */ |
63 |
|
|
*_nbg, /* given background image (optional) */ |
64 |
|
|
*nbg, /* resampled background image */ |
65 |
|
|
*nrgbaD; /* rgba input as double */ |
66 |
|
|
const char *me; |
67 |
|
|
char *outS, *errS; |
68 |
|
|
double _gamma, contr, cfp, cpow, back[3], *rgbaD, r, g, b, a; |
69 |
|
|
airArray *mop; |
70 |
|
|
int E; |
71 |
|
|
size_t min[3], max[3], sx, sy, pi; |
72 |
|
|
unsigned char *outUC, *bgUC; |
73 |
|
|
NrrdResampleInfo *rinfo; |
74 |
|
|
|
75 |
|
|
me = argv[0]; |
76 |
|
|
mop = airMopNew(); |
77 |
|
|
hestOptAdd(&hopt, "i", "nin", airTypeOther, 1, 1, &nin, NULL, |
78 |
|
|
"input nrrd to composite", NULL, NULL, nrrdHestNrrd); |
79 |
|
|
hestOptAdd(&hopt, "c", "contrast", airTypeDouble, 1, 1, &contr, "0.0", |
80 |
|
|
"contrast to apply to RGB values, before gamma. \"0.0\" " |
81 |
|
|
"means no change, \"1.0\" means thresholding, \"-1.0\" " |
82 |
|
|
"means a complete washout."); |
83 |
|
|
hestOptAdd(&hopt, "cfp", "fixed point", airTypeDouble, 1, 1, &cfp, "0.5", |
84 |
|
|
"component level that doesn't change with contrast"); |
85 |
|
|
hestOptAdd(&hopt, "g", "gamma", airTypeDouble, 1, 1, &_gamma, "1.0", |
86 |
|
|
"gamma to apply to image data, after contrast"); |
87 |
|
|
hestOptAdd(&hopt, "b", "background", airTypeDouble, 3, 3, back, "0 0 0", |
88 |
|
|
"background color to composite against; white is " |
89 |
|
|
"1 1 1, not 255 255 255."); |
90 |
|
|
hestOptAdd(&hopt, "bi", "nbg", airTypeOther, 1, 1, &_nbg, "", |
91 |
|
|
"8-bit RGB background image to composite against", |
92 |
|
|
NULL, NULL, nrrdHestNrrd); |
93 |
|
|
hestOptAdd(&hopt, "o", "filename", airTypeString, 1, 1, &outS, |
94 |
|
|
NULL, "file to write output PPM image to"); |
95 |
|
|
hestParseOrDie(hopt, argc-1, argv+1, NULL, me, overInfo, |
96 |
|
|
AIR_TRUE, AIR_TRUE, AIR_TRUE); |
97 |
|
|
airMopAdd(mop, hopt, (airMopper)hestOptFree, airMopAlways); |
98 |
|
|
airMopAdd(mop, hopt, (airMopper)hestParseFree, airMopAlways); |
99 |
|
|
|
100 |
|
|
if (!(3 == nin->dim && 4 <= nin->axis[0].size)) { |
101 |
|
|
fprintf(stderr, "%s: doesn't look like an RGBA nrrd\n", me); |
102 |
|
|
airMopError(mop); return 1; |
103 |
|
|
} |
104 |
|
|
if (nrrdTypeBlock == nin->type) { |
105 |
|
|
fprintf(stderr, "%s: can't use a %s nrrd\n", me, |
106 |
|
|
airEnumStr(nrrdType, nrrdTypeBlock)); |
107 |
|
|
airMopError(mop); return 1; |
108 |
|
|
} |
109 |
|
|
|
110 |
|
|
sx = nin->axis[1].size; |
111 |
|
|
sy = nin->axis[2].size; |
112 |
|
|
if (_nbg) { |
113 |
|
|
if (!(3 == _nbg->dim |
114 |
|
|
&& 3 == _nbg->axis[0].size |
115 |
|
|
&& 2 <= _nbg->axis[1].size |
116 |
|
|
&& 2 <= _nbg->axis[2].size |
117 |
|
|
&& nrrdTypeUChar == _nbg->type)) { |
118 |
|
|
fprintf(stderr, "%s: background not an 8-bit RGB image\n", me); |
119 |
|
|
airMopError(mop); return 1; |
120 |
|
|
} |
121 |
|
|
nbg = nrrdNew(); |
122 |
|
|
airMopAdd(mop, nbg, (airMopper)nrrdNuke, airMopAlways); |
123 |
|
|
if (sx == _nbg->axis[1].size && sy == _nbg->axis[2].size) { |
124 |
|
|
/* no resampling needed, just copy */ |
125 |
|
|
E = nrrdCopy(nbg, _nbg); |
126 |
|
|
} else { |
127 |
|
|
/* have to resample background image to fit. */ |
128 |
|
|
/* because we're using the old resampler, we have to kill off any |
129 |
|
|
space direction information, which is incompatible with |
130 |
|
|
setting per-axis min and max, as is required by the old resampler */ |
131 |
|
|
nrrdOrientationReduce(_nbg, _nbg, AIR_FALSE); |
132 |
|
|
rinfo = nrrdResampleInfoNew(); |
133 |
|
|
airMopAdd(mop, rinfo, (airMopper)nrrdResampleInfoNix, airMopAlways); |
134 |
|
|
rinfo->kernel[0] = NULL; |
135 |
|
|
nrrdKernelParse(&(rinfo->kernel[1]), rinfo->parm[1], "tent"); |
136 |
|
|
rinfo->min[1] = _nbg->axis[1].min = 0; |
137 |
|
|
rinfo->max[1] = _nbg->axis[1].max = _nbg->axis[1].size-1; |
138 |
|
|
rinfo->samples[1] = sx; |
139 |
|
|
nrrdKernelParse(&(rinfo->kernel[2]), rinfo->parm[2], "tent"); |
140 |
|
|
rinfo->min[2] = _nbg->axis[2].min = 0; |
141 |
|
|
rinfo->max[2] = _nbg->axis[2].max = _nbg->axis[2].size-1; |
142 |
|
|
rinfo->samples[2] = sy; |
143 |
|
|
rinfo->renormalize = AIR_TRUE; |
144 |
|
|
rinfo->round = AIR_TRUE; |
145 |
|
|
E = nrrdSpatialResample(nbg, _nbg, rinfo); |
146 |
|
|
} |
147 |
|
|
if (E) { |
148 |
|
|
fprintf(stderr, "%s: trouble:\n%s", me, errS = biffGetDone(NRRD)); |
149 |
|
|
free(errS); return 1; |
150 |
|
|
} |
151 |
|
|
} else { |
152 |
|
|
nbg = NULL; |
153 |
|
|
} |
154 |
|
|
|
155 |
|
|
ninD = nrrdNew(); |
156 |
|
|
airMopAdd(mop, ninD, (airMopper)nrrdNuke, airMopAlways); |
157 |
|
|
nrgbaD = nrrdNew(); |
158 |
|
|
airMopAdd(mop, nrgbaD, (airMopper)nrrdNuke, airMopAlways); |
159 |
|
|
nout=nrrdNew(); |
160 |
|
|
airMopAdd(mop, nout, (airMopper)nrrdNuke, airMopAlways); |
161 |
|
|
|
162 |
|
|
E = 0; |
163 |
|
|
if (nrrdTypeIsIntegral[nin->type]) { |
164 |
|
|
if (!E) E |= nrrdUnquantize(ninD, nin, nrrdTypeDouble); |
165 |
|
|
} else if (nrrdTypeFloat == nin->type) { |
166 |
|
|
if (!E) E |= nrrdConvert(ninD, nin, nrrdTypeDouble); |
167 |
|
|
} else { |
168 |
|
|
if (!E) E |= nrrdCopy(ninD, nin); |
169 |
|
|
} |
170 |
|
|
min[0] = min[1] = min[2] = 0; |
171 |
|
|
max[0] = 3; |
172 |
|
|
max[1] = sx-1; |
173 |
|
|
max[2] = sy-1; |
174 |
|
|
if (!E) E |= nrrdCrop(nrgbaD, ninD, min, max); |
175 |
|
|
if (!E) E |= nrrdPPM(nout, sx, sy); |
176 |
|
|
if (E) { |
177 |
|
|
fprintf(stderr, "%s: trouble:\n%s", me, errS = biffGetDone(NRRD)); |
178 |
|
|
free(errS); return 1; |
179 |
|
|
} |
180 |
|
|
|
181 |
|
|
contr = AIR_CLAMP(-1, contr, 1); |
182 |
|
|
cpow = tan(AIR_AFFINE(-1.000001, contr, 1.000001, 0, AIR_PI/2)); |
183 |
|
|
outUC = (unsigned char*)nout->data; |
184 |
|
|
bgUC = nbg ? (unsigned char *)nbg->data : NULL; |
185 |
|
|
rgbaD = (double *)nrgbaD->data; |
186 |
|
|
for (pi=0; pi<sx*sy; pi++) { |
187 |
|
|
r = AIR_CLAMP(0, rgbaD[0], 1); |
188 |
|
|
g = AIR_CLAMP(0, rgbaD[1], 1); |
189 |
|
|
b = AIR_CLAMP(0, rgbaD[2], 1); |
190 |
|
|
a = AIR_CLAMP(0, rgbaD[3], 1); |
191 |
|
|
if (1 != cpow) { |
192 |
|
|
r = docontrast(r, cfp, cpow); |
193 |
|
|
g = docontrast(g, cfp, cpow); |
194 |
|
|
b = docontrast(b, cfp, cpow); |
195 |
|
|
} |
196 |
|
|
r = pow(r, 1.0/_gamma); |
197 |
|
|
g = pow(g, 1.0/_gamma); |
198 |
|
|
b = pow(b, 1.0/_gamma); |
199 |
|
|
if (bgUC) { |
200 |
|
|
r = a*r + (1-a)*bgUC[0 + 3*pi]/255; |
201 |
|
|
g = a*g + (1-a)*bgUC[1 + 3*pi]/255; |
202 |
|
|
b = a*b + (1-a)*bgUC[2 + 3*pi]/255; |
203 |
|
|
} else { |
204 |
|
|
r = a*r + (1-a)*back[0]; |
205 |
|
|
g = a*g + (1-a)*back[1]; |
206 |
|
|
b = a*b + (1-a)*back[2]; |
207 |
|
|
} |
208 |
|
|
outUC[0] = airIndex(0.0, r, 1.0, 256); |
209 |
|
|
outUC[1] = airIndex(0.0, g, 1.0, 256); |
210 |
|
|
outUC[2] = airIndex(0.0, b, 1.0, 256); |
211 |
|
|
outUC += 3; |
212 |
|
|
rgbaD += 4; |
213 |
|
|
} |
214 |
|
|
|
215 |
|
|
if (nrrdSave(outS, nout, NULL)) { |
216 |
|
|
fprintf(stderr, "%s: trouble:\n%s", me, errS = biffGetDone(NRRD)); |
217 |
|
|
free(errS); return 1; |
218 |
|
|
} |
219 |
|
|
|
220 |
|
|
airMopOkay(mop); |
221 |
|
|
return 0; |
222 |
|
|
} |