1 #if 1 /*######################################################################*/
2 #include <stdlib.h>
3 #include <stdio.h>
4
main(int argc,char * argv[])5 int main( int argc , char *argv[] ) /* Obituary: 27 Apr 2021 */
6 {
7 printf("\n"
8 "*******************************************************************\n"
9 "Program 3dNwarpCalc has been retired, and is no longer available :(\n"
10 "*******************************************************************\n"
11 ) ;
12 exit(0) ;
13 }
14
15 #else /*######################################################################*/
16 #include "mrilib.h"
17 #include "r_new_resam_dset.h"
18
19 #ifdef USE_OMP
20 #include <omp.h>
21 #endif
22
23 #include "mri_genalign.c"
24 #include "mri_genalign_util.c"
25 #include "mri_nwarp.c"
26
27 /*----------------------------------------------------------------------------*/
28
NWC_help(void)29 void NWC_help(void)
30 {
31 printf(
32 "Usage: 3dNwarpCalc [options] expression\n"
33 "------\n"
34 " * This program performs calculations on 3D warps defined on a grid.\n"
35 "\n"
36 " * The fundamental idea is a 'stack' of warps, with operators being\n"
37 " applied to the top element(s) of the stack.\n"
38 " ++ If you don't know what a computer science 'stack' is, see\n"
39 " http://en.wikipedia.org/wiki/Stack_(data_structure)\n"
40 " ++ Also see 1dmatcalc for a similar implementation of a stack\n"
41 " of matrix operations.\n"
42 " ++ In the explanations below, the stack will be denoted as\n"
43 " [ A B C ... ]\n"
44 " where A is the top element, B the next element, etc.\n"
45 " Operations take place using the top one or two elements.\n"
46 "\n"
47 " * The expression is a single string enclosed in quotes ' or \",\n"
48 " with operators separated by spaces. See EXAMPLES below.\n"
49 " ++ The expression is the last thing on the command line!\n"
50 " ++ For scripting convenience, you can actually break the\n"
51 " expression into multiple strings (after all the '-' options),\n"
52 " and they will be re-assembled into one big string for\n"
53 " parsing and processing.\n"
54 "\n"
55 "** Note that to get any output, you will have to use the '&write'\n"
56 " operator at least once. Otherwise, the program computes stuff\n"
57 " and then just throws it away. (Fun perhaps, but useless.)\n"
58 "\n"
59 " * To actually use a 3D warp to transform a dataset, you must run the\n"
60 " program 3dNwarpApply:\n"
61 " ++ Note that the warp used in 3dNwarpApply does not need to be on the\n"
62 " same grid as the dataset being transformed. You can define a warp\n"
63 " on a high-resolution anatomical grid and apply it to a low-resolution\n"
64 " functional dataset, for example -- 3dNwarpApply will figure it out.\n"
65 " ++ On the other hand, all the warps in 3dNwarpCalc must be defined on\n"
66 " the same spatial grid!\n"
67 " ++ (LATER) You can use the &apply command to transform a 3D dataset\n"
68 " from within 3dNwarpCalc, if you don't need the special capabilities\n"
69 " of 3dNwarpApply.\n"
70 "\n"
71 " * Operations such as &invert and &sqrt may produce artifacts and be\n"
72 " inaccurate near the edges of the 3D grid, since they might require\n"
73 " extrapolating the warp to the outside of the grid in places, where\n"
74 " there is no information.\n"
75 "\n"
76 " * For convenience, you can break the expression up into multiple strings\n"
77 " (after all the options), and they will be re-assembled into the single\n"
78 " expression string the controlling C function needs to work.\n"
79 " ++ But note that since the '&' and '()' characters are special to the shell\n"
80 " you have to put the expression string(s) inside single quote ' or double\n"
81 " quote \" pairs, or else ugly things will happen.\n"
82 " (Not as ugly as having a hippopotamus step on your head, but almost.)\n"
83 "\n"
84 "OPTIONS\n"
85 "-------\n"
86 " -interp iii == 'iii' is the interpolation mode:\n"
87 " ++ Modes allowed are a subset of those in 3dAllineate:\n"
88 " linear quintic wsinc5\n"
89 " ++ The default interpolation mode is 'quintic'.\n"
90 " ++ 'linear' is much faster but less accurate.\n"
91 " ++ 'wsinc5' is much slower but more accurate.\n"
92 "\n"
93 " -ainterp jjj == 'jjj' is the interpolation mode for the '&apply' operation.\n"
94 " ++ Modes allowed here are\n"
95 " NN linear cubic quintic wsinc5\n"
96 " ++ If this option isn't given, then the value from '-interp'\n"
97 " is used (which should be good enough for government work).\n"
98 "\n"
99 " -verb == print (to stderr) various fun messages along the road\n"
100 " ++ A second '-verb' gives you even more fun!\n"
101 "\n"
102 "BUT WHERE DO WARPS COME FROM, MOMMY?\n"
103 "------------------------------------\n"
104 " * The program 3dAllineate with the -nwarp_save option will save a\n"
105 " displacement representation of a nonlinear warp to a 3D dataset\n"
106 " with 3 sub-bricks (1 for each of x, y, and z).\n"
107 "\n"
108 " * The contents of these sub-bricks are the displacments of each voxel in mm.\n"
109 " ++ The identity warp would be all zero, for example.\n"
110 "\n"
111 " * An input warp dataset can contain extra sub-bricks -- only the first 3\n"
112 " are used.\n"
113 "\n"
114 #if 0
115 " * Warp datasets output by this program have a 4th sub-brick, labeled\n"
116 " 'hexvol', which contains the volume of each distorted hexahedron in\n"
117 " the grid. This sub-brick is NOT used in any application of the\n"
118 " warp, such as 3dNwarpApply or further runs of 3dNwarpCalc, but can\n"
119 " help you understand how much distortion is present in the warp.\n"
120 " ++ If there are any negative values in the hexvol sub-brick, this\n"
121 " indicates that something bad happened in the warp calculation.\n"
122 #else
123 " * If you want the volume distortion at each voxel, use the program\n"
124 " 3dNwarpFuncs.\n"
125 #endif
126 "\n"
127 "OPERATORS\n"
128 "---------\n"
129 " * In the explanations below, the single character 'x' represents a 3D\n"
130 " coordinate vector, and a capital letter such as 'A' represents a\n"
131 " whole 3D warp function, whose output at a particular location is 'A(x)'.\n"
132 " * You can replace the '&' character that starts a command with '%%' or '@',\n"
133 " if that is more convenient for you.\n"
134 " * Operator names are not case sensitive: &INVERT is the same as &invert.\n"
135 "\n"
136 "&readnwarp(FF) == Read a 3D warp from a file and place it on top of the stack.\n"
137 " *OR* The input file should be a 3D dataset with 3 sub-bricks\n"
138 "&readwarp(FF) (volumes) storing the xyz displacments of each grid point.\n"
139 "\n"
140 "&identwarp(FF) == Create an identity warp (all displacements 0) on the grid\n"
141 " of a 3D dataset specified by the filename 'FF'.\n"
142 " ++ This operation is to be used to create a starting point\n"
143 " for calculations that otherwise do not involve a warp\n"
144 " defined on a grid, such a polynomial warps.\n"
145 " ++ The actual data in 'FF' is ignored by '&identwarp'; only the\n"
146 " 3D grid definition in the header is actually needed.\n"
147 " ----**==>> ++ Either '&identwarp' or '&readnwarp' should be the first\n"
148 " operation, in order to define the grid for all subsequent\n"
149 " calculations.\n"
150 "\n"
151 "&readpoly(FF) == The input is a text file with one line of numbers\n"
152 " specifying a warp as a polynomial, as output from\n"
153 " '3dAllineate -1Dparam_save'.\n"
154 " ++ The count of values determines the type of warp:\n"
155 " 12 ==> affine (shifts+angles+scales+shears)\n"
156 " 64 ==> cubic (3rd order) polyomial\n"
157 " 172 ==> quintic (5th order) polynomial\n"
158 " 364 ==> heptic (7th order) polynomial\n"
159 " 664 ==> nonic (9th order) polynomial\n"
160 " ++ Any other count of values on the single input line is\n"
161 " illegal, unconstitutional, against the laws of God,\n"
162 " fattening, and will make you get red pimples on your nose.\n"
163 " ++ The parameters could come, most probably, from using\n"
164 " 3dAllineate with the '-1Dparam_save' and '-nwarp' options.\n"
165 "\n"
166 "&read4x4(FF) == Read an affine 4x4 transform matrix directly; the input\n"
167 " file should contain 12 numbers in the order:\n"
168 " r11 r12 r13 r14 r21 r22 r23 r24 r31 r32 r33 r34\n"
169 " which will be organized into the 3D transformation matrix:\n"
170 " r11 r12 r13 r14\n"
171 " r21 r22 r23 r24\n"
172 " r31 r32 r33 r34\n"
173 " 0.0 0.0 0.0 1.0\n"
174 " ++ This matrix defines the transformation from input spatial\n"
175 " DICOM coordinates (x,y,z) to output coordinates, in mm.\n"
176 " ++ One way to get this matrix is via '3dAllineate -1Dmatrix_save'.\n"
177 " ++ This matrix should have non-zero determinant!\n"
178 "\n"
179 "&write(FF) == Write the 3D warp on the top of the stack to a file.\n"
180 " The output file is always in a 3D nwarp (dataset) format\n"
181 " -- NEVER a matrix or polynomial.\n"
182 "\n"
183 "&dup == Push the duplicate of the top of the stack onto the stack:\n"
184 " [ A B C ... ] goes to [ A A B C ... ] after &dup.\n"
185 "\n"
186 "&swap == Interchange the top two elements of the stack:\n"
187 " [ A B C ... ] goes to [ B A C ... ] after &swap\n"
188 " ++ You can swap other elements of the stack by using\n"
189 " indexes in the form '&swap(p,q)' where 'p' and 'q'\n"
190 " are distinct non-negative integers indicating depth\n"
191 " into the stack; '&swap' is equivalent to '&swap(0,1)'.\n"
192 "\n"
193 "&pop == Remove (and delete) the top element from the stack:\n"
194 " [ A B C ... ] goes to [ B C ... ] after &pop\n"
195 "\n"
196 "&compose == If the stack is [ A(x) B(x) C(x) ... ], compute the warp\n"
197 " *OR* B(A(x)) and replace these top 2 elements with the result:\n"
198 "&mult [ A B C ... ] goes to [ B(A(x)) C(x) ... ] after &compose\n"
199 " ++ If you wanted to compute A(B(x)), then you would use the\n"
200 " operator combination '&swap &compose'.\n"
201 "\n"
202 "&invert == Replace top element of the stack with its inverse:\n"
203 " the warp J(x) such that A(J(x)) = x.\n"
204 " ++ Inversion is done via a functional iteration:\n"
205 " Jnew(x) = Jold( 2*x - A(Jold(x)) )\n"
206 " which requires 1 warp composition and 1 warp interpolation\n"
207 " for each step.\n"
208 " ++ &invert and &invsqrt (and thus &sqrt) are slow operations\n"
209 " due to the iterative nature of the calculations.\n"
210 #ifdef USE_OMP
211 " ++ Multiple CPUS (via OpenMP) are used to help speed up\n"
212 " these functions.\n"
213 #else
214 " ++ On a system with OpenMP enabled, multiple CPUs would\n"
215 " be used to speed up these functions. However, this\n"
216 " binary copy of 3dNwarpCalc has not been compiled with\n"
217 " OpenMP (alas).\n"
218 #endif
219 " ++ The '-verb' option to 3dNwarpCalc will show you the\n"
220 " progress of the iterations for &invert and &invsqrt.\n"
221 "\n"
222 "&sqrt == Replace top element of the stack with its 'square root':\n"
223 " the warp Q(x) such that Q(Q(x)) = A(x).\n"
224 " ++ NOTE: not all warps have square roots, so this operation\n"
225 " is not guaranteed to work. Be careful out there.\n"
226 " ++ Nor is the square root of a nonlinear operator guaranteed\n"
227 " to be unique!\n"
228 #ifndef USE_SQRTPAIR
229 " ++ The basic algorithm used computes the inverse of Q(x),\n"
230 " as in &invsqrt, so an extra warp inversion is required\n"
231 " at the end here, so this operation is slower than &invsqrt.\n"
232 #endif
233 "\n"
234 "&invsqrt == Replace the top element of the stack with the inverse of\n"
235 " its square root: the warp R(x) such that A(R(R(x)) = x.\n"
236 " ++ '&sqrtinv' is a synonym for this operation, since I always\n"
237 " have trouble remembering which one is correct-imundo-ific.\n"
238 #ifndef USE_SQRTPAIR
239 " ++ This operation is based on the functional iteration\n"
240 " Rnew(x) = Rold( 1.5*x - 0.5*A(Rold(Rold(x))) )\n"
241 " which is adapted from the Schulz iteration for matrix\n"
242 " square roots, and requires 2 warp compositions and 1 warp\n"
243 " interpolation for each step.\n"
244 #else
245 " ++ This operation is based on a functional iteration\n"
246 " adapted from the Denman-Beavers method for computing\n"
247 " the square root of a matrix:\n"
248 " initialize Y(x) = A(x) and Z(x) = x; then iterate\n"
249 " Ynew(x) = 0.5*(Yold(x)+inv(Zold(x)))\n"
250 " Znew(x) = 0.5*(Zold(x)+inv(Yold(x)))\n"
251 " which converges to Y=sqrt(A) and Z=invsqrt(A).\n"
252 " ++ For speed, these square root iterations are always done\n"
253 " with linear interpolation, no matter what '-interp' is.\n"
254 #endif
255 "\n"
256 "&sqrtpair == Compute both &sqrtinv and &sqrt, and leave both of them\n"
257 " on the stack -- &sqrt on top, &sqrtinv 'below' it.\n"
258 "\n"
259 "&sqr == Replace the top element of the stack with its 'square':\n"
260 " the warp S(x) = A(A(x)). Equivalent to '&dup &compose'.\n"
261 " ++ To compute the fourth power of a warp: '&sqr &sqr'\n"
262 " ++ To compute the third power of a warp: '&dup &sqr &compose'\n"
263 " ++ '&square' is a synonym for this operation.\n"
264 "\n"
265 "&scale(a) == Scale the top-of-stack warp displacements by numerical\n"
266 " factor 'a' in all 3 dimensions.\n"
267 " ++ NOTE: this might make the warp non-invertible, (e.g., give\n"
268 " negative results in 'hexvol') for large enough 'a'.\n"
269 " Proceed at your own risk!\n"
270 " ++ If a=0, then the result is the identity warp, since\n"
271 " all the displacements are now 0.\n"
272 " ++ The case a=-1 is NOT the inverse warp!\n"
273 "\n"
274 "&sum == Add the displacements of the two warps on the stack,\n"
275 " then replace BOTH of them with the result.\n"
276 " ++ You can do something like '&sum(0.5,0.5)' to average\n"
277 " the displacements, or '&sum(1,-1)' to difference them.\n"
278 " In this case, the first value scales the displacements\n"
279 " of the stack's top warp, and the second value scales the\n"
280 " displacements of the stack's second warp.\n"
281 " ++ NOTE: you can produce a non-invertible warp this way!\n"
282 "\n"
283 "&apply(DD,PP) == Apply the 3D warp at the top of the stack to a dataset\n"
284 " whose name is given by the 'DD' argument, to produce\n"
285 " a dataset whose prefix is given by the 'PP' argument.\n"
286 " ++ This operation does not affect the stack of warps.\n"
287 " ++ &apply is provided to make your life simpler and happier :-)\n"
288 " ++ &apply is like 3dNwarpApply with the output dataset PP\n"
289 " always being on the same grid as the input dataset DD.\n"
290 " ++ The grid of dataset DD does NOT have to be the same as the\n"
291 " grid defining the warp defined on the stack. If needed,\n"
292 " the warp will be interpolated to be used with DD.\n"
293 " ++ Program 3dNwarpApply provides more options to control\n"
294 " the way that a warp is applied to a dataset; for example,\n"
295 " to control the output grid spacing.\n"
296 "\n"
297 "EXAMPLES\n"
298 "--------\n"
299 "** Read a warp from a dataset, invert it, save the inverse.\n"
300 "\n"
301 " 3dNwarpCalc '&readnwarp(Warp+tlrc.HEAD) &invert &write(WarpInv)'\n"
302 "\n"
303 "** Do the same, but also compute the composition of the warp with the inverse,\n"
304 " and save that -- ideally, the output warp displacements would be identically\n"
305 " zero (i.e., the identity warp), and the 'hexvol' entries would be constant\n"
306 " and equal to the voxel volume.\n"
307 "\n"
308 " 3dNwarpCalc -verb '&readnwarp(Warp+tlrc.HEAD) &dup &invert' \\\n"
309 " '&write(WarpInv) &compose &write(WarpOut)'\n"
310 "\n"
311 "** Read in a warp, compute its inverse square root, then square that, and compose\n"
312 " the result with the original warp -- the result should be the identity warp\n"
313 " (i.e., all zero displacments) -- except for numerical errors, of course.\n"
314 "\n"
315 " 3dNwarpCalc '&readnwarp(Warp+tlrc.HEAD) &dup &invsqrt &sqr &compose &write(WarpOut)'\n"
316 ) ;
317
318 printf(
319 "\n"
320 "AUTHOR -- RWCox -- August 2011\n"
321 ) ;
322
323 PRINT_AFNI_OMP_USAGE("3dNwarpCalc",NULL) ; PRINT_COMPILE_DATE ;
324 exit(0) ;
325 }
326
327 /*----------------------------------------------------------------------------*/
328 /* This program is basically a wrapper for function NwarpCalcRPN() */
329 /*----------------------------------------------------------------------------*/
330
main(int argc,char * argv[])331 int main( int argc , char *argv[] )
332 {
333 int iarg=1 , verb=0 , interp_code=MRI_QUINTIC , ainterp_code=-666 ;
334 char *expr ; int nexpr , narg ;
335 THD_3dim_dataset *oset ;
336
337 AFNI_SETUP_OMP(0) ; /* 24 Jun 2013 */
338
339 if( argc < 2 || strcasecmp(argv[1],"-help") == 0 ) NWC_help() ;
340
341 mainENTRY("3dNwarpCalc"); machdep();
342 AFNI_logger("3dNwarpCalc",argc,argv);
343 PRINT_VERSION("3dNwarpCalc"); AUTHOR("Bob the Warped");
344 (void)COX_clock_time() ;
345
346 while( iarg < argc && argv[iarg][0] == '-' ){
347
348 if( strcasecmp(argv[iarg],"-NN") == 0 || strncasecmp(argv[iarg],"-nearest",6) == 0 ){
349 WARNING_message("NN interpolation not legal here -- switched to linear") ;
350 interp_code = MRI_LINEAR ; iarg++ ; continue ;
351 }
352 if( strncasecmp(argv[iarg],"-linear",4)==0 || strncasecmp(argv[iarg],"-trilinear",6)==0 ){
353 interp_code = MRI_LINEAR ; iarg++ ; continue ;
354 }
355 if( strncasecmp(argv[iarg],"-cubic",4)==0 || strncasecmp(argv[iarg],"-tricubic",6)==0 ){
356 WARNING_message("cubic interplation not legal here -- switched to quintic") ;
357 interp_code = MRI_QUINTIC ; iarg++ ; continue ;
358 }
359 if( strncasecmp(argv[iarg],"-quintic",4)==0 || strncasecmp(argv[iarg],"-triquintic",6)==0 ){
360 interp_code = MRI_QUINTIC ; iarg++ ; continue ;
361 }
362 if( strncasecmp(argv[iarg],"-wsinc",5) == 0 ){
363 interp_code = MRI_WSINC5 ; iarg++ ; continue ;
364 }
365 if( strncasecmp(argv[iarg],"-interp",5)==0 ){
366 char *inam ;
367 if( ++iarg >= argc ) ERROR_exit("no argument after '%s' :-(",argv[iarg-1]) ;
368 inam = argv[iarg] ; if( *inam == '-' ) inam++ ;
369 if( strcasecmp(inam,"NN")==0 || strncasecmp(inam,"nearest",5)==0 ){
370 WARNING_message("NN interpolation not legal here -- changed to linear") ;
371 interp_code = MRI_LINEAR ;
372 } else if( strncasecmp(inam,"linear",3)==0 || strncasecmp(inam,"trilinear",5)==0 ){
373 interp_code = MRI_LINEAR ;
374 } else if( strncasecmp(inam,"cubic",3)==0 || strncasecmp(inam,"tricubic",5)==0 ){
375 WARNING_message("cubic interplation not legal here -- changed to quintic") ;
376 interp_code = MRI_QUINTIC ;
377 } else if( strncasecmp(inam,"quintic",3)==0 || strncasecmp(inam,"triquintic",5)==0 ){
378 interp_code = MRI_QUINTIC ;
379 } else if( strncasecmp(inam,"wsinc",4)==0 ){
380 interp_code = MRI_WSINC5 ;
381 } else {
382 ERROR_exit("Unknown code '%s' after '%s' :-(",argv[iarg],argv[iarg-1]) ;
383 }
384 iarg++ ; continue ;
385 }
386
387 /*---------------*/
388
389 if( strncasecmp(argv[iarg],"-ainterp",6)==0 ){
390 char *inam ;
391 if( ++iarg >= argc ) ERROR_exit("no argument after '%s' :-(",argv[iarg-1]) ;
392 inam = argv[iarg] ; if( *inam == '-' ) inam++ ;
393 if( strcasecmp(inam,"NN")==0 || strncasecmp(inam,"nearest",5)==0 ){
394 ainterp_code = MRI_NN ;
395 } else if( strncasecmp(inam,"linear",3)==0 || strncasecmp(inam,"trilinear",5)==0 ){
396 ainterp_code = MRI_LINEAR ;
397 } else if( strncasecmp(inam,"cubic",3)==0 || strncasecmp(inam,"tricubic",5)==0 ){
398 ainterp_code = MRI_CUBIC ;
399 } else if( strncasecmp(inam,"quintic",3)==0 || strncasecmp(inam,"triquintic",5)==0 ){
400 ainterp_code = MRI_QUINTIC ;
401 } else if( strncasecmp(inam,"wsinc",4)==0 ){
402 ainterp_code = MRI_WSINC5 ;
403 } else {
404 ERROR_exit("Unknown code '%s' after '%s' :-(",argv[iarg],argv[iarg-1]) ;
405 }
406 iarg++ ; continue ;
407 }
408
409 /*---------------*/
410
411 if( strcasecmp(argv[iarg],"-verb") == 0 ){
412 verb++ ; iarg++ ; continue ;
413 }
414
415 /*---------------*/
416
417 ERROR_message("Bizarre and Unknown option '%s' :-(",argv[iarg]) ;
418 suggest_best_prog_option(argv[0],argv[iarg]) ;
419 exit(1) ;
420
421 }
422
423 if( iarg >= argc ) ERROR_exit("No command line expression :-(") ;
424
425 /*--- Assemble all remaining args into the expression ---*/
426
427 expr = strdup(argv[iarg++]) ;
428 for( ; iarg < argc ; iarg++ ){
429 nexpr = strlen(expr) ;
430 narg = strlen(argv[iarg]) ; if( narg == 0 ) continue ;
431 expr = (char *)realloc( expr , sizeof(char)*(nexpr+narg+4) ) ;
432 strcat(expr," ") ; strcat(expr,argv[iarg]) ;
433 }
434
435 if( ainterp_code < 0 ) ainterp_code = interp_code ;
436
437 /*--- All the work is done herein ---*/
438
439 NwarpCalcRPN_verb(verb) ;
440
441 oset = NwarpCalcRPN( expr , NULL , interp_code , ainterp_code ) ;
442
443 /*--- run away screaming into the night, never to be seen again ---*/
444
445 free(expr) ; if( oset != NULL ) DSET_delete(oset) ;
446
447 INFO_message("total CPU time = %.1f sec Elapsed = %.1f\n",
448 COX_cpu_time() , COX_clock_time() ) ;
449
450 exit(0) ;
451 }
452 #endif /*#####################################################################*/
453