Home My Page Projects Code Snippets Project Openings diderot
Summary Activity Tracker Tasks SCM

SCM Repository

[diderot] Annotation of /branches/lamont/src/dnorm/dnorm.c
ViewVC logotype

Annotation of /branches/lamont/src/dnorm/dnorm.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2081 - (view) (download) (as text)

1 : jhr 114 /*
2 :     Teem: Tools to process and visualize scientific data and images
3 :     Copyright (C) 2008, 2007, 2006, 2005 Gordon Kindlmann
4 :     Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998 University of Utah
5 :    
6 :     This library is free software; you can redistribute it and/or
7 :     modify it under the terms of the GNU Lesser General Public License
8 :     (LGPL) as published by the Free Software Foundation; either
9 :     version 2.1 of the License, or (at your option) any later version.
10 :     The terms of redistributing and/or modifying this software also
11 :     include exceptions to the LGPL that facilitate static linking.
12 :    
13 :     This library is distributed in the hope that it will be useful,
14 :     but WITHOUT ANY WARRANTY; without even the implied warranty of
15 :     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 :     Lesser General Public License for more details.
17 :    
18 :     You should have received a copy of the GNU Lesser General Public License
19 :     along with this library; if not, write to Free Software Foundation, Inc.,
20 :     51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
21 :     */
22 :    
23 :     #ifdef DIDEROT
24 :     # include <teem/air.h>
25 :     # include <teem/biff.h>
26 :     # include <teem/nrrd.h>
27 :     #else
28 :     # include "../nrrd.h"
29 :     #endif
30 :    
31 :     char *dnormInfo = ("Normalizes nrrd representation for Diderot. "
32 :     "Forces information about kind and orientation into "
33 :     "a consistent form, and nixes various other fields. ");
34 :    
35 :     int
36 : glk 1234 main(int argc, const char **argv) {
37 : glk 1235 const char *me;
38 :     char *outS;
39 : jhr 114 hestOpt *hopt;
40 :     hestParm *hparm;
41 :     airArray *mop;
42 :    
43 :     char *err;
44 :     Nrrd *nin, *nout;
45 :     NrrdIoState *nio;
46 : glk 1877 int kindIn, kindOut, headerOnly, haveMM, trivialOrient, recenter;
47 : jhr 114 unsigned int kindAxis, axi, si, sj;
48 : glk 1877 double sscl;
49 : jhr 114
50 :     me = argv[0];
51 :     mop = airMopNew();
52 :     hparm = hestParmNew();
53 :     hopt = NULL;
54 :     airMopAdd(mop, hparm, (airMopper)hestParmFree, airMopAlways);
55 : glk 125 hestOptAdd(&hopt, "h,header", NULL, airTypeInt, 0, 0, &headerOnly, NULL,
56 :     "output header of nrrd file only, not the data itself");
57 : glk 1739 hestOptAdd(&hopt, "to", NULL, airTypeInt, 0, 0, &trivialOrient, NULL,
58 :     "(*t*rivial *o*rientation) "
59 :     "even if the input nrrd comes with full orientation or "
60 :     "per-axis min-max info, ignore it and instead assert the "
61 :     "most trivial mapping between index and world space");
62 : glk 1877 hestOptAdd(&hopt, "c", NULL, airTypeInt, 0, 0, &recenter, NULL,
63 :     "re-locate output spaceOrigin so that field is centered "
64 :     "around origin of space coordinates");
65 :     hestOptAdd(&hopt, "s", "scl", airTypeDouble, 1, 1, &sscl, "1.0",
66 :     "when contriving orientation information, distance between "
67 :     "samples to use");
68 : jhr 114 hestOptAdd(&hopt, "i", "nin", airTypeOther, 1, 1, &nin, NULL,
69 :     "input image", NULL, NULL, nrrdHestNrrd);
70 :     hestOptAdd(&hopt, "o", "nout", airTypeString, 1, 1, &outS, "-",
71 :     "output filename", NULL);
72 :    
73 :     hestParseOrDie(hopt, argc-1, argv+1, hparm,
74 :     me, dnormInfo, AIR_TRUE, AIR_TRUE, AIR_TRUE);
75 :     airMopAdd(mop, hopt, (airMopper)hestOptFree, airMopAlways);
76 :     airMopAdd(mop, hopt, (airMopper)hestParseFree, airMopAlways);
77 :    
78 :     /* can't deal with block type */
79 :     if (nrrdTypeBlock == nin->type) {
80 :     fprintf(stderr, "%s: can only have scalar kinds (not %s)\n", me,
81 :     airEnumStr(nrrdType, nrrdTypeBlock));
82 :     airMopError(mop); exit(1);
83 :     }
84 :    
85 :     /* make sure all kinds are set to something */
86 :     /* see if there's a range kind, verify that there's only one */
87 : jhr 1116 /* set haveMM */
88 :     haveMM = AIR_TRUE;
89 : jhr 114 kindIn = nrrdKindUnknown;
90 :     kindAxis = 0;
91 :     for (axi=0; axi<nin->dim; axi++) {
92 :     if (nrrdKindUnknown == nin->axis[axi].kind
93 :     || nrrdKindIsDomain(nin->axis[axi].kind)) {
94 : jhr 1116 haveMM &= AIR_EXISTS(nin->axis[axi].min);
95 :     haveMM &= AIR_EXISTS(nin->axis[axi].max);
96 :     } else {
97 :     if (nrrdKindUnknown != kindIn) {
98 :     fprintf(stderr, "%s: got non-domain kind %s on axis %u, but already "
99 :     "have %s from axis %u\n", me,
100 :     airEnumStr(nrrdKind, nin->axis[axi].kind), axi,
101 :     airEnumStr(nrrdKind, kindIn), kindAxis);
102 :     airMopError(mop); exit(1);
103 :     }
104 :     kindIn = nin->axis[axi].kind;
105 :     kindAxis = axi;
106 : jhr 114 }
107 :     }
108 :     /* see if the non-domain kind is something we can interpret as a tensor */
109 :     if (nrrdKindUnknown != kindIn) {
110 :     switch (kindIn) {
111 :     /* ======= THESE are the kinds that we can possibly output ======= */
112 :     case nrrdKind2Vector:
113 :     case nrrdKind3Vector:
114 :     case nrrdKind4Vector:
115 :     case nrrdKind2DSymMatrix:
116 :     case nrrdKind2DMatrix:
117 :     case nrrdKind3DSymMatrix:
118 :     case nrrdKind3DMatrix:
119 :     /* =============================================================== */
120 :     kindOut = kindIn;
121 :     break;
122 :     case nrrdKind3Color:
123 :     case nrrdKindRGBColor:
124 :     kindOut = nrrdKind3Vector;
125 :     break;
126 :     case nrrdKind4Color:
127 :     case nrrdKindRGBAColor:
128 :     kindOut = nrrdKind4Vector;
129 :     break;
130 :     default:
131 :     fprintf(stderr, "%s: got non-conforming kind %s on axis %u\n", me,
132 :     airEnumStr(nrrdKind, kindIn), kindAxis);
133 :     airMopError(mop); exit(1);
134 :     break;
135 :     }
136 :     } else {
137 :     kindOut = nrrdKindUnknown;
138 :     }
139 :    
140 : glk 147 /* initialize output by copying */
141 : jhr 114 nout = nrrdNew();
142 :     airMopAdd(mop, nout, (airMopper)nrrdNuke, airMopAlways);
143 :     if (nrrdCopy(nout, nin)) {
144 :     airMopAdd(mop, err = biffGet(NRRD), airFree, airMopAlways);
145 :     fprintf(stderr, "%s: trouble copying:\n%s", me, err);
146 :     airMopError(mop); exit(1);
147 :     }
148 :    
149 :     /* no comments, either advertising the format URL or anything else */
150 :     nio = nrrdIoStateNew();
151 :     airMopAdd(mop, nio, (airMopper)nrrdIoStateNix, airMopAlways);
152 :     nio->skipFormatURL = AIR_TRUE;
153 : glk 125 if (headerOnly) {
154 :     nio->skipData = AIR_TRUE;
155 :     }
156 : jhr 114 nrrdCommentClear(nout);
157 :    
158 :     /* no measurement frame */
159 :     for (si=0; si<NRRD_SPACE_DIM_MAX; si++) {
160 :     for (sj=0; sj<NRRD_SPACE_DIM_MAX; sj++) {
161 :     nout->measurementFrame[si][sj] = AIR_NAN;
162 :     }
163 :     }
164 :    
165 :     /* no key/value pairs */
166 :     nrrdKeyValueClear(nout);
167 :    
168 :     /* no content field */
169 :     nout->content = airFree(nout->content);
170 :    
171 :     /* normalize domain kinds to "space" */
172 :     /* turn off centers (perhaps Diderot should assume cell-centered) */
173 :     /* turn off thickness */
174 :     /* turn off labels and units */
175 :     for (axi=0; axi<nout->dim; axi++) {
176 :     if (nrrdKindUnknown == kindOut) {
177 :     nout->axis[axi].kind = nrrdKindSpace;
178 :     } else {
179 :     nout->axis[axi].kind = (kindAxis == axi
180 :     ? kindOut
181 :     : nrrdKindSpace);
182 :     }
183 :     nout->axis[axi].center = nrrdCenterUnknown;
184 :     nout->axis[axi].thickness = AIR_NAN;
185 :     nout->axis[axi].label = airFree(nout->axis[axi].label);
186 :     nout->axis[axi].units = airFree(nout->axis[axi].units);
187 :     nout->axis[axi].min = AIR_NAN;
188 :     nout->axis[axi].max = AIR_NAN;
189 :     nout->axis[axi].spacing = AIR_NAN;
190 :     }
191 :    
192 :     /* logic of orientation definition:
193 :     if space dimension is known:
194 :     set origin to zero if not already set
195 :     set space direction to unit vector if not already set
196 :     else:
197 :     set origin to zero and all space directions to units
198 :     might be nice to use gage's logic for mapping from world to index,
199 :     but we have to accept a greater variety of kinds and dimensions
200 :     than gage ever has to process.
201 :     */
202 : glk 1739 if (nout->spaceDim && !trivialOrient) {
203 : jhr 114 int saxi = 0;
204 : glk 147 /* we use only the space dimension, not any named space */
205 :     nout->space = nrrdSpaceUnknown;
206 : jhr 114 if (!nrrdSpaceVecExists(nout->spaceDim, nout->spaceOrigin)) {
207 :     nrrdSpaceVecSetZero(nout->spaceOrigin);
208 :     }
209 :     for (axi=0; axi<nout->dim; axi++) {
210 :     if (nrrdKindUnknown == kindOut || kindAxis != axi) {
211 :     if (!nrrdSpaceVecExists(nout->spaceDim,
212 :     nout->axis[axi].spaceDirection)) {
213 :     nrrdSpaceVecSetZero(nout->axis[axi].spaceDirection);
214 : glk 1877 nout->axis[axi].spaceDirection[saxi] = sscl;
215 : jhr 114 }
216 :     saxi++;
217 :     } else {
218 :     nrrdSpaceVecSetNaN(nout->axis[axi].spaceDirection);
219 :     }
220 :     }
221 : glk 1739 } else if (haveMM && !trivialOrient) {
222 : jhr 1116 int saxi = 0;
223 :     for (axi=0; axi<nout->dim; axi++) {
224 :     if (nrrdKindUnknown == kindOut || kindAxis != axi) {
225 :     nrrdSpaceVecSetZero(nout->axis[axi].spaceDirection);
226 :     nout->axis[axi].spaceDirection[saxi]
227 :     = (nin->axis[axi].max - nin->axis[axi].min)/(nin->axis[axi].size-1);
228 :     nout->spaceOrigin[saxi] = nin->axis[axi].min;
229 :     saxi++;
230 :     } else {
231 :     nrrdSpaceVecSetNaN(nout->axis[axi].spaceDirection);
232 :     }
233 :     }
234 :     nout->spaceDim = saxi;
235 : jhr 114 } else {
236 : glk 1740 /* either trivialOrient, or not spaceDim, or not not haveMM */
237 : glk 1743 int saxi = 0;
238 : glk 1740 nout->space = nrrdSpaceUnknown;
239 : jhr 114 nrrdSpaceVecSetZero(nout->spaceOrigin);
240 :     for (axi=0; axi<nout->dim; axi++) {
241 :     if (nrrdKindUnknown == kindOut || kindAxis != axi) {
242 :     nrrdSpaceVecSetZero(nout->axis[axi].spaceDirection);
243 : glk 125 nout->axis[axi].spaceDirection[saxi]
244 :     = (AIR_EXISTS(nin->axis[axi].spacing)
245 :     ? nin->axis[axi].spacing
246 : glk 1877 : sscl);
247 : jhr 114 saxi++;
248 :     } else {
249 :     nrrdSpaceVecSetNaN(nout->axis[axi].spaceDirection);
250 :     }
251 :     }
252 :     nout->spaceDim = saxi;
253 :     }
254 : jhr 1116
255 :     /* space dimension has to match the number of domain axes */
256 : jhr 114 if (nout->dim != nout->spaceDim + !!kindOut) {
257 :     fprintf(stderr, "%s: output dim %d != spaceDim %d + %d %s%s%s\n",
258 :     me, nout->dim, nout->spaceDim, !!kindOut,
259 :     kindOut ? "for non-scalar (" : "(scalar data)",
260 :     kindOut ? airEnumStr(nrrdKind, kindOut) : "",
261 :     kindOut ? ") data" : "");
262 :     airMopError(mop); exit(1);
263 :     }
264 :    
265 : glk 1877 if (recenter) {
266 :     /* sets field's origin so field is centered on the origin. capiche? */
267 :     /* this code was tacked on later than the stuff above, so its
268 :     logic could probably be moved up there, but it seems cleaner to
269 :     have it as a separate post-process */
270 :     double mean[NRRD_SPACE_DIM_MAX];
271 :     nrrdSpaceVecSetZero(mean);
272 :     for (axi=0; axi<nout->dim; axi++) {
273 :     if (nrrdKindUnknown == kindOut || kindAxis != axi) {
274 :     nrrdSpaceVecScaleAdd2(mean, 1.0, mean,
275 :     0.5*(nout->axis[axi].size - 1),
276 :     nout->axis[axi].spaceDirection);
277 :     }
278 :     }
279 :     nrrdSpaceVecScaleAdd2(mean, 1.0, mean,
280 :     1.0, nout->spaceOrigin);
281 :     /* now mean is the center of the field */
282 :     nrrdSpaceVecScaleAdd2(nout->spaceOrigin,
283 :     1.0, nout->spaceOrigin,
284 :     -1.0, mean);
285 :     }
286 :    
287 : jhr 114 if (nrrdSave(outS, nout, nio)) {
288 :     airMopAdd(mop, err = biffGet(NRRD), airFree, airMopAlways);
289 :     fprintf(stderr, "%s: trouble saving \"%s\":\n%s",
290 :     me, outS, err);
291 :     airMopError(mop); exit(1);
292 :     }
293 :    
294 :     airMopOkay(mop);
295 :     exit(0);
296 :     }

root@smlnj-gforge.cs.uchicago.edu
ViewVC Help
Powered by ViewVC 1.0.0