*** src/brain/smooth_missing.c Wed Mar 24 17:56:51 1999 --- /SLOTH-U2/welling/fiasco/development/Fiasco4.2/src/brain//smooth_missing.c Wed Mar 3 18:56:40 1999 *************** *** 0 **** --- 1,207 ---- + /************************************************************ + * * + * smooth_missing.c * + * * + * Permission is hereby granted to any individual or * + * institution for use, copying, or redistribution of * + * this code and associated documentation, provided * + * that such code and documentation are not sold for * + * profit and the following copyright notice is retained * + * in the code and documentation: * + * Copyright (c) 1999 Department of Statistics, * + * Carnegie Mellon University * + * * + * Original programming by Joel Welling 3/1999 * + ************************************************************/ + + #include + #include + #include + #include + #include "mri.h" + #include "fmri.h" + #include "stdcrg.h" + #include "misc.h" + + static char rcsid[] = "$Id: smooth_missing.c,v 1.1 1999/03/03 23:56:23 welling Exp $"; + + + static void smooth_ts( float* corr_ts, float* ts, unsigned char** missing, + int z, int dt ) + { + int t; + int t_start; + int t_end; + int loop; + int in_run= 0; + + for (t=0; t 1 ) && !strcmp( argv[1], "-help" ) ) + { + if( argc == 2 ) + Help( "selecttopic" ); + else + Help( argv[2] ); + } + + + /*** Parse command line ***/ + + cl_scan( argc, argv ); + + /* Get filenames */ + cl_get( "dataout|d", "%option %s[%]", ".dat", outfile ); + cl_get( "headerout|h", "%option %s[%]", "smooth_missing.mri", hdrfile ); + cl_get( "input|i", "%option %s[%]", "input.mri", infile ); + + if (cl_cleanup_check()) { + int i; + fprintf(stderr,"%s: invalid argument in command line:\n ",argv[0]); + for (i=0; i data.root root + data.path path + + + data.file $$ + data.format pgh + + data.file data + data.format format + data.file + + + # + # Experimental Paradigm + # + + conditions 1 + design.images 1 + design.aligned 1 + fixed 0 + IAI ****IAI**** + + design + # condition stimulus.start stimulus.len + 0 0 ****NUM_IMAGES**** + + + + # + # Fixed Parameters + # + + shape.params 4 + drift.degree 3 + paste.form 0 + paste.scale 0.001 + + drift.degree degree + + # + # Prior Specification + # + + prior.baseline 2000.0 1.0e-6 + prior.drift 0.001 1.0 5.0 + prior.resp 0.005 50.0 + prior.sigmasq 1.8 200.0 + + prior.shape + 2.0 4.30 + 4.0 1.0 + 2.0 4.30 + 4.0 0.43 + + + prior.drift.scaled false + + + # + # Algorithm Parameters --- vmpfx + # + + vmpfx.iterations 2500 + vmpfx.evaluations 25000 + vmpfx.scd.tol 1.0e-5 + vmpfx.cov.method analytic + vmpfx.knots ****NUM_KNOTS**** + vmpfx.fix.knots true + + vmpfx.init auto + vmpfx.produce residuals + vmpfx.model full + + + # + # Algorithm Parameters --- vfpfx + # + + vfpfx.samples 10500 + vfpfx.burnin 500 + vfpfx.increment 1 + + vfpfx.init auto + + vfpfx.jump + # name default [value] + baseline 0.5 + shape 0.05 + responsiveness 0.001 + + + vfpfx.se.expand 2.5 + + vfpfx.knots.max 6 + vfpfx.knots.init $ + vfpfx.fix.knots $ + + vfpfx.products + # name class which + baseline cumulants 2 + responsiveness samples 1 + drift samples 25 + shape quantiles 5 + noise.precision quantiles 5 + + + + # + # Output Format + # + + output.root vmpfx_results + output.path "" + output.root oroot + output.path opath + + output.file $$.out + + output.append $<0> + vmpfx false + vfpfx true + true + + + output.file output + output.logfile logfile + *** src/csh/spline_detrend.csh Wed Mar 24 17:24:04 1999 --- /SLOTH-U2/welling/fiasco/development/Fiasco4.2/src/csh//spline_detrend.csh Wed Mar 24 16:32:31 1999 *************** *** 0 **** --- 1,113 ---- + #!/bin/csh -efx + # spline_detrend.csh + #/************************************************************ + # * * + # * Permission is hereby granted to any individual or * + # * institution for use, copying, or redistribution of * + # * this code and associated documentation, provided * + # * that such code and documentation are not sold for * + # * profit and the following copyright notice is retained * + # * in the code and documentation: * + # * Copyright (c) 1999 Department of Statistics, * + # * Carnegie Mellon University * + # * * + # * Original programming by Joel Welling * + # ************************************************************/ + # + echo '#'`date` $0 + echo '#$Id: spline_detrend.csh,v 1.2 1999/03/24 21:32:31 welling Exp $' + if(! -d stat)mkdir stat + if(! -d ps) mkdir ps + if (`printmrifield -i $1 -f $F_DETREND_CHUNK.dimensions` == "xyzt") then + remap -file $1 -chunk $F_DETREND_CHUNK -order vxyzt + endif + permute -memlimit 32000000 -input $1.mri -headerout ${1}_p.mri \ + -dataout .dat -chunk $F_DETREND_CHUNK -order vtxyz + echo '#'`date` + smooth_missing -i ${1}_p.mri -h ${1}_s.mri + rm ${1}_p.mri ${1}_p.dat + echo '#'`date` + + # Generate vmpfx input parameter file for null model + generate_vmpfx_in -proto $F_VMPFX_NPROTO -nimage $F_NIMAGE -nslice \ + $F_NSLICE -droot ${1:t}_s -dpath ${1:h}/ \ + -iai $F_IAI -null > $F_VMPFX_NINPUT + + # Run vmpfx + echo '#'`date` + vmpfx $F_VMPFX_NINPUT + echo '#'`date` + rm ${1}_s.mri ${1}_s.dat + + # Move vmpfx results to appropriate directories + mv vmpfx_results.out out/spline_detrend_vmpfx.out + mv vmpfx_results.log out/spline_detrend_vmpfx.log + mv vmpfx_results.bin ${3}_d.dat + + # Reality check the vmpfx output + set ofile = out/spline_detrend_vmpfx.out + set tot_size = `awk '/max.posterior/ { print $2 ; exit }' < $ofile` + set res_units = `awk '/record/ { print $2 ; exit }' < $ofile` + set res_count = `awk '/residuals/ { print $2 ; exit }' < $ofile` + set res_size = `awk '/residuals/ { print $3 ; exit }' < $ofile` + @ prod = $res_count * $res_size * $res_units + if ( $res_size != 8 ) then + echo 'Unexpected data type for vmpfx output!' + exit -1 + endif + if ( $tot_size != $prod ) then + echo 'Unexpected format for vmpfx output!' + exit -1 + endif + + # Prepare a .mri file for the output + set images_ext_x = `printmrifield -input $1 -field images.extent.x` + set images_ext_y = `printmrifield -input $1 -field images.extent.y` + set images_ext_z = `printmrifield -input $1 -field images.extent.z` + set images_ext_v = `printmrifield -input $1 -field images.extent.v` + set images_ext_t = `printmrifield -input $1 -field images.extent.t` + set images_dims = `printmrifield -input $1 -field images.dimensions` + echo '\!format = pgh' > ${3}_d.mri + echo '\!version = 1.0' >> ${3}_d.mri + echo 'images = [chunk]' >> ${3}_d.mri + echo 'images.datatype = float64' >> ${3}_d.mri + echo 'images.dimensions=' $images_dims >> ${3}_d.mri + echo 'images.file= .dat' >> ${3}_d.mri + echo 'images.extent.x=' $images_ext_x >> ${3}_d.mri + echo 'images.extent.y=' $images_ext_y >> ${3}_d.mri + echo 'images.extent.z=' $images_ext_z >> ${3}_d.mri + echo 'images.extent.t=' $images_ext_t >> ${3}_d.mri + echo 'images.extent.v= 1' >> ${3}_d.mri + echo 'images.offset= 0' >> ${3}_d.mri + echo 'images.size=' `(ls -l ${3}_d.dat | awk '{ print $5 }')` >> ${3}_d.mri + + # Restore missing data + mri_copy_chunk -chunk missing -input $1 -output ${3}_d + + # Convert the output to type float, and throw out the double version + mri_type_convert -float ${3}_d ${3}_f + rm ${3}_d.mri ${3}_d.dat + + # Come up with a mean dataset, and add it back into the residuals + # It's actually faster to do the permutation than to interpolate + # in time series order. + mri_subsample -d t -e 1 -mean $1 ${3}_m + mri_remap -file ${3}_m -order vxyzt + mri_interp -c -d t -e $images_ext_t ${3}_m ${3}_l + rm ${3}_m.mri ${3}_m.dat + permute -memlimit 32000000 -input ${3}_l -headerout ${3}_M -order vtxyz + rm ${3}_l.mri ${3}_l.dat + mri_rpn_math -o $3 '$1,$2,+' ${3}_f ${3}_M + rm ${3}_f.mri ${3}_f.dat ${3}_M.mri ${3}_M.dat + + # Finally done with input; clean up + rm $1.mri $1.dat + + # Permute the output back to image order + echo '#'`date` + permute -memlimit 32000000 -input $3.mri -headerout $2.mri \ + -dataout .dat -chunk $F_DETREND_CHUNK -order vxyzt + echo '#'`date` + + + *** src/fmri/chirprot.c Wed Mar 24 17:24:04 1999 --- /SLOTH-U2/welling/fiasco/development/Fiasco4.2/src/fmri//chirprot.c Thu Oct 29 18:11:43 1998 *************** *** 0 **** --- 1,250 ---- + /************************************************************ + * * + * fshrot.c * + * * + * Permission is hereby granted to any individual or * + * institution for use, copying, or redistribution of * + * this code and associated documentation, provided * + * that such code and documentation are not sold for * + * profit and the following copyright notice is retained * + * in the code and documentation: * + * * + * Original programming by Mark Fitzgerald 2-95 * + * 5-96: Pittsburgh Format, Mark Fitzgerald * + * 10-98: Chirp transformation, Bill Eddy * + ************************************************************/ + + /************************************************************* + * fourier_shift_rot performs a translation and rotation * + * of a complex image using the Fourier shift theorem * + * and the chirp transform for rotations. + * For details, see "Improved Registration using * + * Fourier Interpolation" by Eddy, Fitzgerald, Noll, * + * Journal of Magnetic Resonance in Medicine (1996) * + * and "Two and Three Dimensional Image Rotation Using * + * the FFT" by Cox and Tong (Note that there are five * + * typos on page 3 of Cox and Tong including a critical * + * one in the definition of Z) * + * par is the RegPars structure containing the amount * + * of translation and rotation to apply * + * orig_image is the 2D, complex array containing the * + * image to be translated and rotated * + * moved_image is the returned 2D, complex array containing * + * translated, rotated version of orig_image * + * nx and ny and the slowest-varying and fastest-varying * + * dimensions of orig_image (and moved_image) respectively * + * domain indicates the domain of the image to which the * + * movement is to be applied: * + * 'f' or 'k' if orig_image is a Fourier (k-space) image * + * 't' or 'i' if orig_image is a regular (i-space) image * + *************************************************************/ + + #include + #include + #include + #include + #include "mri.h" + #include "fmri.h" + #include "stdcrg.h" + + static char rcsid[] = "$Id: chirprot.c,v 1.1 1998/10/29 23:11:23 welling Exp $"; + + void fourier_shift_rot( RegPars par, FComplex** orig_image, + FComplex** moved_image, long nx, long ny, + char domain ) + { + short fourier_flag; + FComplex **temp_image1 = NULL, **temp_image2 = NULL, **temp_image3 = NULL; + FComplex **temp_image4; + long x, y, halfx, halfy, ls, mr, p, q; + float angle, temp_angle, cos_ang, sin_ang, alpha, beta; + + /* Determine domain */ + switch( domain ) + { + case 'i': case 'I': case 't': case 'T': + fourier_flag = 0; + break; + case 'f': case 'F': case 'k': case 'K': + fourier_flag = 1; + break; + default: + Error( "Domain argument to fourier_shift_rot unrecognized: (%c).", + domain ); + } + + /* + Message( "#Domain %c, \n", domain ); + */ + + /* Calculate midpoints */ + halfx = nx / 2; + halfy = ny / 2; + + /* Allocate space for intermediate images */ + temp_image1 = Matrix( ny, nx, FComplex ); + temp_image2 = Matrix( ny, nx, FComplex ); + temp_image3 = Matrix( ny, nx, FComplex ); + temp_image4 = Matrix( ny, nx, FComplex ); + + /* Copy original image into local storage */ + memcpy( *temp_image1, *orig_image, + (size_t) ( sizeof(FComplex) * nx * ny ) ); + + if( !fourier_flag ) + { + /* Perform full 2D FFT */ + fft2d( temp_image1, nx, ny, 1, '2', 0, 0 ); + } + + if( par.x_shift || par.y_shift){ + + /* + Message( "#Shift X %f, Y %f \n", par.x_shift, par.y_shift ); + */ + + /* Apply linear phase change to transformed data */ + for( y = 0; y < ny; y++ ) + { + + /* Calculate y-translation contribution to phase change */ + temp_angle = par.y_shift * (float) ( y - halfy ) / (float) ny; + + for( x = 0; x < nx; x++ ) + { + /* Calculate phase shift to apply for particular location */ + angle = - 2.0 * M_PI * (temp_angle + par.x_shift * + (float) ( x - halfx ) / (float) nx); + cos_ang = cos( angle ); + sin_ang = sin( angle ); + + /* Apply phase shift */ + temp_image2[y][x].real = cos_ang * temp_image1[y][x].real - + sin_ang * temp_image1[y][x].imag; + temp_image2[y][x].imag = cos_ang * temp_image1[y][x].imag + + sin_ang * temp_image1[y][x].real; + } + } + } + else { + memcpy( *temp_image2, *temp_image1, + (size_t) ( sizeof(FComplex) * nx * ny ) ); + } + + if( par.rotation ){ + + /* Apply rotation if non-zero */ + + /* + Message( "#Rotation %f \n", par.rotation ); + */ + + /* Chirp transform so we need two inverse DFT's and one */ + /* forward DFT plus a couple of multiplications */ + + /* This trick requires symmetric matrices */ + if (nx != ny) Abort("chirprot error: input matrics not symmetric!"); + + /* Convert rotation from degrees to radians */ + angle = par.rotation * M_PI / 180.0; + + /* Find alpha and beta*/ + alpha = cos( angle )/(float)( nx ); + beta = sin( angle )/(float)( nx ); + + /* + Message( "#Angle %f, Alpha %f, Beta %f \n", angle, alpha, beta ); + */ + + /* First, multiply by Z(ls,mr) */ + for( y = 0; y < ny; y++){ + for( x = 0; x < nx; x++){ + ls = (x - halfx); + mr = (y - halfy); + angle = M_PI * (((float)(ls * ls - mr * mr)) * alpha + + ((float)(2 * ls * mr)) * beta); + cos_ang = cos(angle); + sin_ang = sin(angle); + temp_image1[y][x].real = cos_ang * temp_image2[y][x].real - + sin_ang * temp_image2[y][x].imag; + temp_image1[y][x].imag = cos_ang * temp_image2[y][x].imag + + sin_ang * temp_image2[y][x].real; + + /* Create Z(ls-p,mr+q)* !!Note complex conjugate*/ + temp_image4[y][x].real = cos_ang; + temp_image4[y][x].imag = - sin_ang; + } + } + + /* + Message( "#Multiply by Z(ls,mr) complete. %d %d \n", x, y ); + */ + + /* Convolve these two images */ + + /* First, perform full 2D inverse FFT */ + fft2d( temp_image4, nx, ny, -1, '2', 0, 0 ); + + /* Second, perform full 2D inverse FFT */ + fft2d( temp_image1, nx, ny, -1, '2', 0, 0 ); + + /* Second, multiply images together */ + for( y = 0; y < ny; y++){ + for( x = 0; x < nx; x++){ + temp_image3[y][x].real = temp_image1[y][x].real + * temp_image4[y][x].real - temp_image1[y][x].imag + * temp_image4[y][x].imag; + temp_image3[y][x].imag = temp_image1[y][x].imag + * temp_image4[y][x].real + temp_image1[y][x].real + * temp_image4[y][x].imag; + } + } + + /* + Message( "#Multiply for convolution complete. %d %d \n", x, y ); + */ + + /* Third, perform full 2D forward FFT */ + fft2d( temp_image3, nx, ny, 1, '2', 0, 0 ); + + /* Third, multiply by Z(q,p)* !!Note complex conjugate*/ + for( y = 0; y < ny; y++){ + for( x = 0; x < nx; x++){ + q = (y - halfy); + p = (x - halfx); + angle = M_PI * (((float)(q * q - p * p)) * alpha + + ((float)(2 * q * p)) * beta); + cos_ang = cos( angle ); + sin_ang = - sin( angle ); + temp_image2[y][x].real = cos_ang * temp_image3[y][x].real + - sin_ang * temp_image3[y][x].imag; + temp_image2[y][x].imag = cos_ang * temp_image3[y][x].imag + + sin_ang * temp_image3[y][x].real; + } + } + /* + Message( "#Multiply by Z(q,p)* complete. %d %d \n", x, y ); + */ + } + else { + /* Get back from K space */ + fft2d( temp_image2, nx, ny, 1, '2', 0, 0 ); + } + + /* Copy local storage back to original image */ + memcpy( *moved_image, *temp_image2, + (size_t) ( sizeof(FComplex) * nx * ny ) ); + + /* Free local storage */ + FreeMatrix( temp_image1 ); + FreeMatrix( temp_image2 ); + FreeMatrix( temp_image3 ); + FreeMatrix( temp_image4 ); + + /* + Message( "#Done with chirp.\n"); + */ + + return; + + } *** src/fmri/fft.c Wed Mar 24 17:24:04 1999 --- /SLOTH-U2/welling/fiasco/development/Fiasco4.2/src/fmri//fft.c Wed Sep 9 13:50:00 1998 *************** *** 0 **** --- 1,855 ---- + #include + #include + + static char rcsid[] = "$Id: fft.c,v 1.2 1997/11/26 00:36:50 welling Exp $"; + + int cfft (int n, float * c, float * wsave, int isign); + + static int cfft1_ (int * n, float * c, float * ch, float * wa, + int * ifac, int * isign); + static int cffti1_(int * n, float * wa, int * ifac); + + + /* from fftpkg package in cmlib and netlib -- translated by f2c */ + /* and modified by Audris Mockus and Mark Fitzgerald */ + + /* + Public Routine + */ + + /* + * cfft computes the forward or backward complex discrete fourier transform. + * + * Input Parameters: + * + * n The length of the complex sequence c. The method is + * more efficient when n is the product of small primes. + * + * c A complex array of length n which contains the sequence + * + * wsave a real work array which must be dimensioned at least 4n+15 + * in the program that calls cfft. + * isign 1 for transform, -1 for inverse transform. + * A call of cfft with isign = 1 followed by a call of cfft with + * isign = -1 will multiply the sequence by n. + * + * Output Parameters: + * + * c For j=1,...,n + * + * c(j)=the sum from k=1,...,n of + * + * c(k)*exp(-isign*i*(j-1)*(k-1)*2*pi/n) + * + * where i=sqrt(-1) + */ + /* Example how to use + float x1 [] = {9, 0, -1, 0, 1, 0, -5, 0}; + float x [] = {1, 0, 2, 0, -1, 0, 12, 0}; + float ws [4*8+15]; + + main () + { + int i=4; + cfft (4, x1, ws, 1); + for (i = 0; i < 4; i++) + printf("%f %f\n", x1 [i*2], x1 [i*2+1]); + cfft (4, x, ws, 1); + for (i = 0; i < 4; i++) + printf("%f %f\n", x [i*2], x1 [i*2+1]); + cfft (4, x1, ws, -1); + for (i = 0; i < 4; i++) + printf("%f %f\n", x1 [i*2], x1 [i*2+1]); + cfft (4, x, ws, -1); + for (i = 0; i < 4; i++) + printf("%f %f\n", x [i*2], x1 [i*2+1]); + } + */ + + + + int cfft (int n, float * c, float * wsave, int isign) + { + int iw1, iw2; + + /* Parameter adjustments */ + --c; + --wsave; + + /* Function Body */ + if (n != 1) { + iw1 = n + n + 1; + iw2 = iw1 + n + n; + cffti1_(&n, &wsave[iw1], (int *)&wsave[iw2]); + cfft1_(&n, &c[1], &wsave[1], &wsave[iw1], (int *)&wsave[iw2], (int*)&isign); + } + return 0; + } + + + static int cffti1_(int * n, float * wa, int * ifac) + { + /* Initialized data */ + static int ntryh[4] = { 3,4,2,5 }; + + /* System generated locals */ + int i_1, i_2, i_3; + + /* Local variables */ + float argh; + int idot, ntry, i, j; + float argld; + int i1, k1, l1, l2, ib; + float fi; + int ld, ii, nf, ip, nl, nq, nr; + float arg; + int ido, ipm; + float tpi; + + /* Parameter adjustments */ + --wa; + --ifac; + + /* Function Body */ + nl = *n; + nf = 0; + j = 0; + + L101: + ++j; + if (j - 4 <= 0) ntry = ntryh[j - 1]; + else ntry += 2; + L104: + nq = nl / ntry; + nr = nl - ntry * nq; + if (nr != 0) goto L101; + ++nf; + ifac[nf + 2] = ntry; + nl = nq; + if (ntry == 2 && nf != 1) { + i_1 = nf; + for (i = 2; i <= i_1; ++i) { + ib = nf - i + 2; + ifac[ib + 2] = ifac[ib + 1]; + } + ifac[3] = 2; + } + if (nl != 1) goto L104; + + ifac[1] = *n; + ifac[2] = nf; + tpi = 6.28318530717959; + argh = tpi / (float) (*n); + i = 2; + l1 = 1; + i_1 = nf; + for (k1 = 1; k1 <= i_1; ++k1) { + ip = ifac[k1 + 2]; + ld = 0; + l2 = l1 * ip; + ido = *n / l2; + idot = ido + ido + 2; + ipm = ip - 1; + i_2 = ipm; + for (j = 1; j <= i_2; ++j) { + i1 = i; + wa[i - 1] = 1.0; + wa[i] = 0.0; + ld += l1; + fi = 0.0; + argld = (float) ld * argh; + i_3 = idot; + for (ii = 4; ii <= i_3; ii += 2) { + i += 2; + fi += 1.0; + arg = fi * argld; + wa[i - 1] = cos(arg); + wa[i] = sin(arg); + } + if (ip > 5) { + wa[i1 - 1] = wa[i - 1]; + wa[i1] = wa[i]; + } + } + l1 = l2; + } + return 0; + } + + static int pass_(nac, ido, ip, l1, idl1, cc, c1, c2, ch, ch2, wa, isign) + int *nac; + int *ido, *ip, *l1, *idl1; + float *cc, *c1, *c2, *ch, *ch2, *wa; + int *isign; + { + /* System generated locals */ + int ch_dim1, ch_dim2, ch_offset, cc_dim1, cc_dim2, cc_offset, c1_dim1, + c1_dim2, c1_offset, c2_dim1, c2_offset, ch2_dim1, ch2_offset, + i_1, i_2, i_3; + + /* Local variables */ + int idij, idlj, idot, ipph, i, j, k, l, jc, lc, ik, idj, idl, inc, idp; + float wai, war; + int ipp2; + + /* Parameter adjustments */ + cc_dim1 = *ido; + cc_dim2 = *ip; + cc_offset = cc_dim1 * (cc_dim2 + 1) + 1; + cc -= cc_offset; + c1_dim1 = *ido; + c1_dim2 = *l1; + c1_offset = c1_dim1 * (c1_dim2 + 1) + 1; + c1 -= c1_offset; + c2_dim1 = *idl1; + c2_offset = c2_dim1 + 1; + c2 -= c2_offset; + ch_dim1 = *ido; + ch_dim2 = *l1; + ch_offset = ch_dim1 * (ch_dim2 + 1) + 1; + ch -= ch_offset; + ch2_dim1 = *idl1; + ch2_offset = ch2_dim1 + 1; + ch2 -= ch2_offset; + --wa; + + /* Function Body */ + idot = *ido / 2; + ipp2 = *ip + 2; + ipph = (*ip + 1) / 2; + idp = *ip * *ido; + + if (*ido >= *l1) { + i_1 = ipph; + for (j = 2; j <= i_1; ++j) { + jc = ipp2 - j; + i_2 = *l1; + for (k = 1; k <= i_2; ++k) { + i_3 = *ido; + for (i = 1; i <= i_3; ++i) { + ch[i + (k + j * ch_dim2) * ch_dim1] = + cc[i + (j + k * cc_dim2) * cc_dim1] + + cc[i + (jc + k * cc_dim2) * cc_dim1]; + ch[i + (k + jc * ch_dim2) * ch_dim1] = + cc[i + (j + k * cc_dim2) * cc_dim1] + - cc[i + (jc + k * cc_dim2) * cc_dim1]; + } + } + } + i_1 = *l1; + for (k = 1; k <= i_1; ++k) { + i_2 = *ido; + for (i = 1; i <= i_2; ++i) { + ch[i + (k + ch_dim2) * ch_dim1] = cc[i + (k * cc_dim2 + 1) * cc_dim1]; + } + } + } + else { + i_1 = ipph; + for (j = 2; j <= i_1; ++j) { + jc = ipp2 - j; + i_2 = *ido; + for (i = 1; i <= i_2; ++i) { + i_3 = *l1; + for (k = 1; k <= i_3; ++k) { + ch[i + (k + j * ch_dim2) * ch_dim1] = + cc[i + (j + k * cc_dim2) * cc_dim1] + + cc[i + (jc + k * cc_dim2) * cc_dim1]; + ch[i + (k + jc * ch_dim2) * ch_dim1] = + cc[i + (j + k * cc_dim2) * cc_dim1] + - cc[i + (jc + k * cc_dim2) * cc_dim1]; + } + } + } + i_1 = *ido; + for (i = 1; i <= i_1; ++i) { + i_2 = *l1; + for (k = 1; k <= i_2; ++k) { + ch[i + (k + ch_dim2) * ch_dim1] = cc[i + (k * cc_dim2 + 1) * cc_dim1]; + } + } + } + idl = 2 - *ido; + inc = 0; + i_1 = ipph; + for (l = 2; l <= i_1; ++l) { + lc = ipp2 - l; + idl += *ido; + i_2 = *idl1; + for (ik = 1; ik <= i_2; ++ik) { + c2[ik + l * c2_dim1] = ch2[ik + ch2_dim1] + wa[idl - 1] + * ch2[ik + (ch2_dim1 << 1)]; + c2[ik + lc * c2_dim1] = - *isign * wa[idl] * ch2[ik + *ip * ch2_dim1]; + } + idlj = idl; + inc += *ido; + i_2 = ipph; + for (j = 3; j <= i_2; ++j) { + jc = ipp2 - j; + idlj += inc; + if (idlj > idp) { + idlj -= idp; + } + war = wa[idlj - 1]; + wai = wa[idlj]; + i_3 = *idl1; + for (ik = 1; ik <= i_3; ++ik) { + c2[ik + l * c2_dim1] += war * ch2[ik + j * ch2_dim1]; + c2[ik + lc * c2_dim1] -= *isign * wai * ch2[ik + jc * ch2_dim1]; + } + } + } + i_1 = ipph; + for (j = 2; j <= i_1; ++j) { + i_2 = *idl1; + for (ik = 1; ik <= i_2; ++ik) { + ch2[ik + ch2_dim1] += ch2[ik + j * ch2_dim1]; + } + } + i_1 = ipph; + for (j = 2; j <= i_1; ++j) { + jc = ipp2 - j; + i_2 = *idl1; + for (ik = 2; ik <= i_2; ik += 2) { + ch2[ik - 1 + j * ch2_dim1] = c2[ik - 1 + j * c2_dim1] + - c2[ik + jc * c2_dim1]; + ch2[ik - 1 + jc * ch2_dim1] = c2[ik - 1 + j * c2_dim1] + + c2[ik + jc * c2_dim1]; + ch2[ik + j * ch2_dim1] = c2[ik + j * c2_dim1] + + c2[ik - 1 + jc * c2_dim1]; + ch2[ik + jc * ch2_dim1] = c2[ik + j * c2_dim1] + - c2[ik - 1 + jc * c2_dim1]; + } + } + *nac = 1; + if (*ido != 2) { + *nac = 0; + i_1 = *idl1; + for (ik = 1; ik <= i_1; ++ik) { + c2[ik + c2_dim1] = ch2[ik + ch2_dim1]; + } + i_1 = *ip; + for (j = 2; j <= i_1; ++j) { + i_2 = *l1; + for (k = 1; k <= i_2; ++k) { + c1[(k + j * c1_dim2) * c1_dim1 + 1] = + ch[(k + j * ch_dim2) * ch_dim1 + 1]; + c1[(k + j * c1_dim2) * c1_dim1 + 2] = + ch[(k + j * ch_dim2) * ch_dim1 + 2]; + } + } + if (idot <= *l1) { + idij = 0; + i_1 = *ip; + for (j = 2; j <= i_1; ++j) { + idij += 2; + i_2 = *ido; + for (i = 4; i <= i_2; i += 2) { + idij += 2; + i_3 = *l1; + for (k = 1; k <= i_3; ++k) { + c1[i - 1 + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1] + * ch[i - 1 + (k + j * ch_dim2) * ch_dim1] + *isign * wa[idij] + * ch[i + (k + j * ch_dim2) * ch_dim1]; + c1[i + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1] + * ch[i + (k + j * ch_dim2) * ch_dim1] - *isign * wa[idij] + * ch[i - 1 + (k + j * ch_dim2) * ch_dim1]; + } + } + } + } + else { + idj = 2 - *ido; + i_1 = *ip; + for (j = 2; j <= i_1; ++j) { + idj += *ido; + i_2 = *l1; + for (k = 1; k <= i_2; ++k) { + idij = idj; + i_3 = *ido; + for (i = 4; i <= i_3; i += 2) { + idij += 2; + c1[i - 1 + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1] + * ch[i - 1 + (k + j * ch_dim2) * ch_dim1] + + *isign * wa[idij] * ch[i + (k + j * ch_dim2) * ch_dim1]; + c1[i + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1] + * ch[i + (k + j * ch_dim2) * ch_dim1] - *isign * wa[idij] + * ch[i - 1 + (k + j * ch_dim2) * ch_dim1]; + } + } + } + } + } + return 0; + } + + static int pass2_(ido, l1, cc, ch, wa1, isign) + int *ido, *l1; + float *cc, *ch, *wa1; + int *isign; + { + /* System generated locals */ + int cc_dim1, cc_offset, ch_dim1, ch_dim2, ch_offset, i_1, i_2; + + /* Local variables */ + int i, k; + float ti2, tr2; + + /* Parameter adjustments */ + cc_dim1 = *ido; + cc_offset = cc_dim1 * 3 + 1; + cc -= cc_offset; + ch_dim1 = *ido; + ch_dim2 = *l1; + ch_offset = ch_dim1 * (ch_dim2 + 1) + 1; + ch -= ch_offset; + --wa1; + + /* Function Body */ + if (*ido <= 2) { + i_1 = *l1; + for (k = 1; k <= i_1; ++k) { + ch[(k + ch_dim2) * ch_dim1 + 1] = cc[((k << 1) + 1) * cc_dim1 + 1] + + cc[((k << 1) + 2) * cc_dim1 + 1]; + ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1] = cc[((k << 1) + 1) * cc_dim1 + 1] + - cc[((k << 1) + 2) * cc_dim1 + 1]; + ch[(k + ch_dim2) * ch_dim1 + 2] = cc[((k << 1) + 1) * cc_dim1 + 2] + + cc[((k << 1) + 2) * cc_dim1 + 2]; + ch[(k + (ch_dim2 << 1)) * ch_dim1 + 2] = cc[((k << 1) + 1) * cc_dim1 + 2] + - cc[((k << 1) + 2) * cc_dim1 + 2]; + } + } + else { + i_1 = *l1; + for (k = 1; k <= i_1; ++k) { + i_2 = *ido; + for (i = 2; i <= i_2; i += 2) { + ch[i - 1 + (k + ch_dim2) * ch_dim1] + = cc[i - 1 + ((k << 1) + 1) * cc_dim1] + + cc[i - 1 + ((k << 1) + 2) * cc_dim1]; + tr2 = cc[i - 1 + ((k << 1) + 1) * cc_dim1] + - cc[i - 1 + ((k << 1) + 2) * cc_dim1]; + ch[i + (k + ch_dim2) * ch_dim1] = cc[i + ((k << 1) + 1) * cc_dim1] + + cc[i + ((k << 1) + 2) * cc_dim1]; + ti2 = cc[i + ((k << 1) + 1) * cc_dim1] + - cc[i + ((k << 1) + 2) * cc_dim1]; + ch[i + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * ti2 + - *isign * wa1[i] * tr2; + ch[i - 1 + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * tr2 + + *isign * wa1[i] * ti2; + } + } + } + return 0; + } + + static int pass3_(ido, l1, cc, ch, wa1, wa2, isign) + int *ido, *l1; + float *cc, *ch, *wa1, *wa2; + int *isign; + { + /* System generated locals */ + int cc_dim1, cc_offset, ch_dim1, ch_dim2, ch_offset, i_1, i_2; + + /* Local variables */ + float taui, taur; + int i, k; + float ci2, ci3, di2, di3, cr2, cr3, dr2, dr3, ti2, tr2; + + /* Parameter adjustments */ + cc_dim1 = *ido; + cc_offset = (cc_dim1 << 2) + 1; + cc -= cc_offset; + ch_dim1 = *ido; + ch_dim2 = *l1; + ch_offset = ch_dim1 * (ch_dim2 + 1) + 1; + ch -= ch_offset; + --wa1; + --wa2; + + /* Function Body */ + taur = -.5; + taui = -(*isign) * .866025403784439; + if (*ido == 2) { + i_1 = *l1; + for (k = 1; k <= i_1; ++k) { + tr2 = cc[(k * 3 + 2) * cc_dim1 + 1] + cc[(k * 3 + 3) * cc_dim1 + 1]; + cr2 = cc[(k * 3 + 1) * cc_dim1 + 1] + taur * tr2; + ch[(k + ch_dim2) * ch_dim1 + 1] = cc[(k * 3 + 1) * cc_dim1 + 1] + tr2; + ti2 = cc[(k * 3 + 2) * cc_dim1 + 2] + cc[(k * 3 + 3) * cc_dim1 + 2]; + ci2 = cc[(k * 3 + 1) * cc_dim1 + 2] + taur * ti2; + ch[(k + ch_dim2) * ch_dim1 + 2] = cc[(k * 3 + 1) * cc_dim1 + 2] + ti2; + cr3 = taui * (cc[(k * 3 + 2) * cc_dim1 + 1] + - cc[(k * 3 + 3) * cc_dim1 + 1]); + ci3 = taui * (cc[(k * 3 + 2) * cc_dim1 + 2] + - cc[(k * 3 + 3) * cc_dim1 + 2]); + ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1] = cr2 - ci3; + ch[(k + ch_dim2 * 3) * ch_dim1 + 1] = cr2 + ci3; + ch[(k + (ch_dim2 << 1)) * ch_dim1 + 2] = ci2 + cr3; + ch[(k + ch_dim2 * 3) * ch_dim1 + 2] = ci2 - cr3; + } + } + else { + i_1 = *l1; + for (k = 1; k <= i_1; ++k) { + i_2 = *ido; + for (i = 2; i <= i_2; i += 2) { + tr2 = cc[i - 1 + (k * 3 + 2) * cc_dim1] + + cc[i - 1 + (k * 3 + 3) * cc_dim1]; + cr2 = cc[i - 1 + (k * 3 + 1) * cc_dim1] + taur * tr2; + ch[i - 1 + (k + ch_dim2) * ch_dim1] = cc[i - 1 + (k * 3 + 1) + * cc_dim1] + tr2; + ti2 = cc[i + (k * 3 + 2) * cc_dim1] + cc[i + (k * 3 + 3) * cc_dim1]; + ci2 = cc[i + (k * 3 + 1) * cc_dim1] + taur * ti2; + ch[i + (k + ch_dim2) * ch_dim1] = cc[i + (k * 3 + 1) * cc_dim1] + ti2; + cr3 = taui * (cc[i - 1 + (k * 3 + 2) * cc_dim1] + - cc[i - 1 + (k * 3 + 3) * cc_dim1]); + ci3 = taui * (cc[i + (k * 3 + 2) * cc_dim1] + - cc[i + (k * 3 + 3) * cc_dim1]); + dr2 = cr2 - ci3; + dr3 = cr2 + ci3; + di2 = ci2 + cr3; + di3 = ci2 - cr3; + ch[i + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * di2 + - *isign * wa1[i] * dr2; + ch[i - 1 + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * dr2 + + *isign * wa1[i] * di2; + ch[i + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * di3 + - *isign * wa2[i] * dr3; + ch[i - 1 + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * dr3 + + *isign * wa2[i] * di3; + } + } + } + return 0; + } + + static int pass4_(ido, l1, cc, ch, wa1, wa2, wa3, isign) + int *ido, *l1; + float *cc, *ch, *wa1, *wa2, *wa3; + int *isign; + { + /* System generated locals */ + int cc_dim1, cc_offset, ch_dim1, ch_dim2, ch_offset, i_1, i_2; + + /* Local variables */ + int i, k; + float ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4; + + /* Parameter adjustments */ + cc_dim1 = *ido; + cc_offset = cc_dim1 * 5 + 1; + cc -= cc_offset; + ch_dim1 = *ido; + ch_dim2 = *l1; + ch_offset = ch_dim1 * (ch_dim2 + 1) + 1; + ch -= ch_offset; + --wa1; + --wa2; + --wa3; + + /* Function Body */ + if (*ido == 2) { + i_1 = *l1; + for (k = 1; k <= i_1; ++k) { + ti1 = cc[((k << 2) + 1) * cc_dim1 + 2] + - cc[((k << 2) + 3) * cc_dim1 + 2]; + ti2 = cc[((k << 2) + 1) * cc_dim1 + 2] + + cc[((k << 2) + 3) * cc_dim1 + 2]; + tr4 = *isign * (cc[((k << 2) + 2) * cc_dim1 + 2] + - cc[((k << 2) + 4) * cc_dim1 + 2]); + ti3 = cc[((k << 2) + 2) * cc_dim1 + 2] + + cc[((k << 2) + 4) * cc_dim1 + 2]; + tr1 = cc[((k << 2) + 1) * cc_dim1 + 1] + - cc[((k << 2) + 3) * cc_dim1 + 1]; + tr2 = cc[((k << 2) + 1) * cc_dim1 + 1] + + cc[((k << 2) + 3) * cc_dim1 + 1]; + ti4 = *isign * (cc[((k << 2) + 4) * cc_dim1 + 1] + - cc[((k << 2) + 2) * cc_dim1 + 1]); + tr3 = cc[((k << 2) + 2) * cc_dim1 + 1] + + cc[((k << 2) + 4) * cc_dim1 + 1]; + ch[(k + ch_dim2) * ch_dim1 + 1] = tr2 + tr3; + ch[(k + ch_dim2 * 3) * ch_dim1 + 1] = tr2 - tr3; + ch[(k + ch_dim2) * ch_dim1 + 2] = ti2 + ti3; + ch[(k + ch_dim2 * 3) * ch_dim1 + 2] = ti2 - ti3; + ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1] = tr1 + tr4; + ch[(k + (ch_dim2 << 2)) * ch_dim1 + 1] = tr1 - tr4; + ch[(k + (ch_dim2 << 1)) * ch_dim1 + 2] = ti1 + ti4; + ch[(k + (ch_dim2 << 2)) * ch_dim1 + 2] = ti1 - ti4; + } + } + else { + i_1 = *l1; + for (k = 1; k <= i_1; ++k) { + i_2 = *ido; + for (i = 2; i <= i_2; i += 2) { + ti1 = cc[i + ((k << 2) + 1) * cc_dim1] + - cc[i + ((k << 2) + 3) * cc_dim1]; + ti2 = cc[i + ((k << 2) + 1) * cc_dim1] + + cc[i + ((k << 2) + 3) * cc_dim1]; + ti3 = cc[i + ((k << 2) + 2) * cc_dim1] + + cc[i + ((k << 2) + 4) * cc_dim1]; + tr4 = *isign * (cc[i + ((k << 2) + 2) * cc_dim1] + - cc[i + ((k << 2) + 4) * cc_dim1]); + tr1 = cc[i - 1 + ((k << 2) + 1) * cc_dim1] + - cc[i - 1 + ((k << 2) + 3) * cc_dim1]; + tr2 = cc[i - 1 + ((k << 2) + 1) * cc_dim1] + + cc[i - 1 + ((k << 2) + 3) * cc_dim1]; + ti4 = *isign * (cc[i - 1 + ((k << 2) + 4) * cc_dim1] + - cc[i - 1 + ((k << 2) + 2) * cc_dim1]); + tr3 = cc[i - 1 + ((k << 2) + 2) * cc_dim1] + + cc[i - 1 + ((k << 2) + 4) * cc_dim1]; + ch[i - 1 + (k + ch_dim2) * ch_dim1] = tr2 + tr3; + cr3 = tr2 - tr3; + ch[i + (k + ch_dim2) * ch_dim1] = ti2 + ti3; + ci3 = ti2 - ti3; + cr2 = tr1 + tr4; + cr4 = tr1 - tr4; + ci2 = ti1 + ti4; + ci4 = ti1 - ti4; + ch[i - 1 + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * cr2 + + *isign * wa1[i] * ci2; + ch[i + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * ci2 + - *isign * wa1[i] * cr2; + ch[i - 1 + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * cr3 + + *isign * wa2[i] * ci3; + ch[i + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * ci3 + - *isign * wa2[i] * cr3; + ch[i - 1 + (k + (ch_dim2 << 2)) * ch_dim1] = wa3[i - 1] * cr4 + + *isign * wa3[i] * ci4; + ch[i + (k + (ch_dim2 << 2)) * ch_dim1] = wa3[i - 1] * ci4 + - *isign * wa3[i] * cr4; + } + } + } + return 0; + } + + static int pass5_(ido, l1, cc, ch, wa1, wa2, wa3, wa4, isign) + int *ido, *l1; + float *cc, *ch, *wa1, *wa2, *wa3, *wa4; + int *isign; + { + /* System generated locals */ + int cc_dim1, cc_offset, ch_dim1, ch_dim2, ch_offset, i_1, i_2; + + /* Local variables */ + int i, k; + float ci2, ci3, ci4, ci5, di3, di4, di5, di2, cr2, cr3, cr5, cr4, ti2, ti3, + ti4, ti5, dr3, dr4, dr5, dr2, tr2, tr3, tr4, tr5, ti11, ti12, tr11, + tr12; + + /* Parameter adjustments */ + cc_dim1 = *ido; + cc_offset = cc_dim1 * 6 + 1; + cc -= cc_offset; + ch_dim1 = *ido; + ch_dim2 = *l1; + ch_offset = ch_dim1 * (ch_dim2 + 1) + 1; + ch -= ch_offset; + --wa1; + --wa2; + --wa3; + --wa4; + + /* Function Body */ + tr11 = .309016994374947; + ti11 = -(*isign) * .951056516295154; + tr12 = -.809016994374947; + ti12 = -(*isign) * .587785252292473; + if (*ido == 2) { + i_1 = *l1; + for (k = 1; k <= i_1; ++k) { + ti5 = cc[(k * 5 + 2) * cc_dim1 + 2] - cc[(k * 5 + 5) * cc_dim1 + 2]; + ti2 = cc[(k * 5 + 2) * cc_dim1 + 2] + cc[(k * 5 + 5) * cc_dim1 + 2]; + ti4 = cc[(k * 5 + 3) * cc_dim1 + 2] - cc[(k * 5 + 4) * cc_dim1 + 2]; + ti3 = cc[(k * 5 + 3) * cc_dim1 + 2] + cc[(k * 5 + 4) * cc_dim1 + 2]; + tr5 = cc[(k * 5 + 2) * cc_dim1 + 1] - cc[(k * 5 + 5) * cc_dim1 + 1]; + tr2 = cc[(k * 5 + 2) * cc_dim1 + 1] + cc[(k * 5 + 5) * cc_dim1 + 1]; + tr4 = cc[(k * 5 + 3) * cc_dim1 + 1] - cc[(k * 5 + 4) * cc_dim1 + 1]; + tr3 = cc[(k * 5 + 3) * cc_dim1 + 1] + cc[(k * 5 + 4) * cc_dim1 + 1]; + ch[(k + ch_dim2) * ch_dim1 + 1] = cc[(k * 5 + 1) * cc_dim1 + 1] + + tr2 + tr3; + ch[(k + ch_dim2) * ch_dim1 + 2] = cc[(k * 5 + 1) * cc_dim1 + 2] + + ti2 + ti3; + cr2 = cc[(k * 5 + 1) * cc_dim1 + 1] + tr11 * tr2 + tr12 * tr3; + ci2 = cc[(k * 5 + 1) * cc_dim1 + 2] + tr11 * ti2 + tr12 * ti3; + cr3 = cc[(k * 5 + 1) * cc_dim1 + 1] + tr12 * tr2 + tr11 * tr3; + ci3 = cc[(k * 5 + 1) * cc_dim1 + 2] + tr12 * ti2 + tr11 * ti3; + cr5 = ti11 * tr5 + ti12 * tr4; + ci5 = ti11 * ti5 + ti12 * ti4; + cr4 = ti12 * tr5 - ti11 * tr4; + ci4 = ti12 * ti5 - ti11 * ti4; + ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1] = cr2 - ci5; + ch[(k + ch_dim2 * 5) * ch_dim1 + 1] = cr2 + ci5; + ch[(k + (ch_dim2 << 1)) * ch_dim1 + 2] = ci2 + cr5; + ch[(k + ch_dim2 * 3) * ch_dim1 + 2] = ci3 + cr4; + ch[(k + ch_dim2 * 3) * ch_dim1 + 1] = cr3 - ci4; + ch[(k + (ch_dim2 << 2)) * ch_dim1 + 1] = cr3 + ci4; + ch[(k + (ch_dim2 << 2)) * ch_dim1 + 2] = ci3 - cr4; + ch[(k + ch_dim2 * 5) * ch_dim1 + 2] = ci2 - cr5; + } + } + else { + i_1 = *l1; + for (k = 1; k <= i_1; ++k) { + i_2 = *ido; + for (i = 2; i <= i_2; i += 2) { + ti5 = cc[i + (k * 5 + 2) * cc_dim1] - cc[i + (k * 5 + 5) * cc_dim1]; + ti2 = cc[i + (k * 5 + 2) * cc_dim1] + cc[i + (k * 5 + 5) * cc_dim1]; + ti4 = cc[i + (k * 5 + 3) * cc_dim1] - cc[i + (k * 5 + 4) * cc_dim1]; + ti3 = cc[i + (k * 5 + 3) * cc_dim1] + cc[i + (k * 5 + 4) * cc_dim1]; + tr5 = cc[i - 1 + (k * 5 + 2) * cc_dim1] + - cc[i - 1 + (k * 5 + 5) * cc_dim1]; + tr2 = cc[i - 1 + (k * 5 + 2) * cc_dim1] + + cc[i - 1 + (k * 5 + 5) * cc_dim1]; + tr4 = cc[i - 1 + (k * 5 + 3) * cc_dim1] + - cc[i - 1 + (k * 5 + 4) * cc_dim1]; + tr3 = cc[i - 1 + (k * 5 + 3) * cc_dim1] + + cc[i - 1 + (k * 5 + 4) * cc_dim1]; + ch[i - 1 + (k + ch_dim2) * ch_dim1] = cc[i - 1 + (k * 5 + 1) + * cc_dim1] + tr2 + tr3; + ch[i + (k + ch_dim2) * ch_dim1] = cc[i + (k * 5 + 1) + * cc_dim1] + ti2 + ti3; + cr2 = cc[i - 1 + (k * 5 + 1) * cc_dim1] + tr11 * tr2 + tr12 * tr3; + ci2 = cc[i + (k * 5 + 1) * cc_dim1] + tr11 * ti2 + tr12 * ti3; + + cr3 = cc[i - 1 + (k * 5 + 1) * cc_dim1] + tr12 * tr2 + tr11 * tr3; + ci3 = cc[i + (k * 5 + 1) * cc_dim1] + tr12 * ti2 + tr11 * ti3; + + cr5 = ti11 * tr5 + ti12 * tr4; + ci5 = ti11 * ti5 + ti12 * ti4; + cr4 = ti12 * tr5 - ti11 * tr4; + ci4 = ti12 * ti5 - ti11 * ti4; + dr3 = cr3 - ci4; + dr4 = cr3 + ci4; + di3 = ci3 + cr4; + di4 = ci3 - cr4; + dr5 = cr2 + ci5; + dr2 = cr2 - ci5; + di5 = ci2 - cr5; + di2 = ci2 + cr5; + ch[i - 1 + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * dr2 + + *isign * wa1[i] * di2; + ch[i + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * di2 + - *isign * wa1[i] * dr2; + ch[i - 1 + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * dr3 + + *isign * wa2[i] * di3; + ch[i + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * di3 + - *isign * wa2[i] * dr3; + ch[i - 1 + (k + (ch_dim2 << 2)) * ch_dim1] = wa3[i - 1] * dr4 + + *isign * wa3[i] * di4; + ch[i + (k + (ch_dim2 << 2)) * ch_dim1] = wa3[i - 1] * di4 + - *isign * wa3[i] * dr4; + ch[i - 1 + (k + ch_dim2 * 5) * ch_dim1] = wa4[i - 1] * dr5 + + *isign * wa4[i] * di5; + ch[i + (k + ch_dim2 * 5) * ch_dim1] = wa4[i - 1] * di5 + - *isign * wa4[i] * dr5; + } + } + } + return 0; + } + + + /* + Internal Routines + */ + + static int cfft1_ (int * n, float * c, float * ch, float * wa, + int * ifac, int * isign) + { + /* System generated locals */ + int i_1; + + /* Local variables */ + int idot; + int i, k1, l1, l2, n2, na, nf, ip, iw, ix2, ix3, ix4, nac, ido, idl1; + + /* Parameter adjustments */ + --c; + --ch; + --wa; + --ifac; + + /* Function Body */ + nf = ifac[2]; + na = 0; + l1 = 1; + iw = 1; + i_1 = nf; + for (k1 = 1; k1 <= i_1; ++k1) { + ip = ifac[k1 + 2]; + l2 = ip * l1; + ido = *n / l2; + idot = ido + ido; + idl1 = idot * l1; + if (ip == 4) { + ix2 = iw + idot; + ix3 = ix2 + idot; + if (na == 0) { + pass4_(&idot, &l1, &c[1], &ch[1], &wa[iw], &wa[ix2], &wa[ix3], isign); + } + else { + pass4_(&idot, &l1, &ch[1], &c[1], &wa[iw], &wa[ix2], &wa[ix3], isign); + } + na = 1 - na; + } + else if (ip == 2) { + if (na == 0) { + pass2_(&idot, &l1, &c[1], &ch[1], &wa[iw], isign); + } + else { + pass2_(&idot, &l1, &ch[1], &c[1], &wa[iw], isign); + } + na = 1 - na; + } + else if (ip == 3) { + ix2 = iw + idot; + if (na == 0) { + pass3_(&idot, &l1, &c[1], &ch[1], &wa[iw], &wa[ix2], isign); + } + else { + pass3_(&idot, &l1, &ch[1], &c[1], &wa[iw], &wa[ix2], isign); + } + na = 1 - na; + } + else if (ip == 5) { + ix2 = iw + idot; + ix3 = ix2 + idot; + ix4 = ix3 + idot; + if (na == 0) { + pass5_(&idot, &l1, &c[1], &ch[1], &wa[iw], &wa[ix2], &wa[ix3], + &wa[ix4], isign); + } + else { + pass5_(&idot, &l1, &ch[1], &c[1], &wa[iw], &wa[ix2], &wa[ix3], + &wa[ix4], isign); + } + na = 1 - na; + } + else { + if (na == 0) { + pass_(&nac, &idot, &ip, &l1, &idl1, &c[1], &c[1], &c[1], &ch[1], + &ch[1], &wa[iw], isign); + } + else { + pass_(&nac, &idot, &ip, &l1, &idl1, &ch[1], &ch[1], &ch[1], &c[1], + &c[1], &wa[iw], isign); + } + if (nac != 0) { + na = 1 - na; + } + } + l1 = l2; + iw += (ip - 1) * idot; + } + if (na != 0) { + n2 = *n + *n; + i_1 = n2; + for (i = 1; i <= i_1; ++i) { + c[i] = ch[i]; + } + } + return 0; + } *** src/fmri/fft3d_help.help Wed Mar 24 17:24:04 1999 --- /SLOTH-U2/welling/fiasco/development/Fiasco4.2/src/fmri//fft3d_help.help Wed Feb 24 19:56:36 1999 *************** *** 0 **** --- 1,25 ---- + *FFTAlgorithm + + The FFTW package (available from http://theory.lcs.mit.edu/~fftw) is + used to implement the fast Fourier transforms used in this routine. + FFTW (Fastest Fourier Transform in the West) is distributed under + the Gnu General Public License. + + *FFTEnvironment + + If the environment variable F_WISDOM_FILE is set, its value will be + used as the basis for the name of a file in which FFTW will store + "wisdom". This is machine-specific information about what precise + FFT algorithm is fastest on the specific hardware. This information + can then be reused, saving a couple of seconds the next time a + program using FFTW runs. + + Because this information is machine specific, it is best to assign + F_WISDOM_FILE a machine-specific value, for example + "in/fftw_wisdom.${PVM_ARCH}" . + + In parallel runs, FFTW will create one wisdom file per processor, + by appending processor numbers to the file name. This is necessary + to avoid errors due to multiple processes reading and writing the same + file. + *** src/fmri/fft3d_tester.c Wed Mar 24 17:24:04 1999 --- /SLOTH-U2/welling/fiasco/development/Fiasco4.2/src/fmri//fft3d_tester.c Thu Oct 1 18:09:13 1998 *************** *** 0 **** --- 1,126 ---- + /************************************************************ + * * + * fft3d.c * + * * + * Permission is hereby granted to any individual or * + * institution for use, copying, or redistribution of * + * this code and associated documentation, provided * + * that such code and documentation are not sold for * + * profit and the following copyright notice is retained * + * in the code and documentation: * + * * + * Original programming by Joel Welling, 9/98 * + ************************************************************/ + + /* This utility reads a file containing Pgh MRI data in "xyz" form, + * and does various FFT things to it to test "fft3d". It then + * writes the resulting data out. + */ + + #include + #include + #include + #include + #include "mri.h" + #include "fmri.h" + + int main(int argc, char* argv[]) + { + MRI_Dataset *Input= NULL; + MRI_Dataset *Output= NULL; + float* obuf; + float* ibuf; + FComplex* data; + int nx; + int ny; + int nz; + int i; + int j; + int k; + + if (argc != 3) { + fprintf(stderr,"Usage: %s in.mri out.mri\n", argv[0]); + fprintf(stderr, + " in.mri must be an xyz Pgh MRI dataset in image space,\n"); + fprintf(stderr, + " containing float data.\n"); + } + else { + Input= mri_open_dataset( argv[1], MRI_READ ); + if (!Input) { + fprintf(stderr,"%s: cannot open %s for reading!\n",argv[0],argv[1]); + exit(-1); + } + if (strcmp(mri_get_string(Input,"images.dimensions"),"xyz")) { + fprintf(stderr,"%s: input data must have dims xyz!\n",argv[0]); + exit(-1); + } + if (strcmp(mri_get_string(Input,"images.datatype"),"float32")) { + fprintf(stderr,"%s: input data must have type float32!\n",argv[0]); + exit(-1); + } + + nx= mri_get_int(Input,"images.extent.x"); + ny= mri_get_int(Input,"images.extent.y"); + nz= mri_get_int(Input,"images.extent.z"); + + Output= mri_copy_dataset( argv[2], Input ); + if (!Output) { + fprintf(stderr,"%s: error creating %s for writing!\n",argv[0],argv[2]); + exit(-1); + } + + if (!(obuf= (float*)malloc(nx*ny*nz*sizeof(float)))) { + fprintf(stderr,"Unable to allocate %d floats!\n",nx*ny*nz); + exit(-1); + } + if (!(data= (FComplex*)malloc(nx*ny*nz*sizeof(FComplex)))) { + fprintf(stderr,"Unable to allocate %d FComplex!\n",nx*ny*nz); + exit(-1); + } + + ibuf= mri_get_chunk(Input,"images",nx*ny*nz,0,MRI_FLOAT); + + for (i=0; i3 + 3->xy,z; 3->yz,x; 3->xz,y + 3->x,y,z + x->x + y->y + z->z + xy->xy + xz->xz + y,z->yz + yz->y,z + yz->yz + Failed tests: + */ + + fft3d(data,nx,ny,nz,+1,"yz"); + fft3d(data,nx,ny,nz,-1,"yz"); + + for (i=0; i + #include + #include + #include "mri.h" + #include "fmri.h" + #include "stdcrg.h" + #include "misc.h" + + /* Notes- + * -Need more insight into what makes bad regions bad + * -Note that the routine expects data presented in z-fastest order! + */ + + static char rcsid[] = "$Id: fshrot3d.c,v 1.5 1999/02/26 23:34:45 welling Exp $"; + + /* How bad does cancellation have to get before we consider a quat bad? */ + #define CANCELLATION_TOL 0.001 + + /* How long to look for a good quat */ + #define REPAIR_QUAT_MAX_STEPS 10 + + /* How big of steps to take looking for a good quat */ + #define REPAIR_QUAT_STEPSIZE 0.01 + + /* Accessor for 3D arrays */ + #define MEM(matrix,nx,ny,nz,x,y,z) matrix[((((x)*ny)+(y))*nz)+(z)] + + typedef struct shear_param_struct { + double a, b, c, d, e, f, g, h; + } ShearParams; + + typedef struct trans_param_struct { + double dx, dy, dz; + } TransParams; + + static int debug= 0; + static int count_shear_x= 0; + static int count_shear_y= 0; + static int count_shear_z= 0; + static int count_set_phase= 0; + static int count_calls= 0; + + void fshrot3d_set_debug( int flag ) + { + debug= flag; + } + + void fshrot3d_clear_shear_counts() + { + count_shear_x= count_shear_y= count_shear_z= count_calls= + count_set_phase= 0; + } + + void fshrot3d_get_shear_counts( int* xcount, int* ycount, int* zcount, + int* phase_count, + double *mean_shears_per_call) + { + *xcount= count_shear_x; + *ycount= count_shear_y; + *zcount= count_shear_z; + *phase_count= count_set_phase; + *mean_shears_per_call= + (double)(count_shear_x + count_shear_y + count_shear_z + count_set_phase) + / (double)count_calls; + } + + static void quat_to_shears( Quat* q, ShearParams* s) + { + if (q->x==0.0 && q->y==0.0 && q->z==0.0) { + /* Special case the identity quaternion */ + s->a= s->b= s->c= s->d= s->e= s->f= s->g= s->h= 0.0; + if (debug) { + fprintf(stderr, + "Result of trivial decomp: a=%f, b=%f, c=%f, d= %f\n", + s->a,s->b,s->c,s->d); + fprintf(stderr, + " e= %f, f=%f, g=%f, h=%f\n", + s->e,s->f,s->g,s->h); + } + } + else { + double X= q->x; + double Y= q->y; + double Z= q->z; + double W= q->w; + + double t1= (2.0*X*Y-2.0*Z*W); + double t2= (-(2.0*X*Z-2.0*Y*W)*t1 + +(1.0-2.0*Y*Y-2.0*Z*Z)*(2.0*Y*Z+2.0*X*W)); + + if (debug) fprintf(stderr,"Decomposing quat (%f %f %f %f)\n",X,Y,Z,W); + + s->a= -(-(1.0-2.0*X*X-2.0*Z*Z)*(1.0-2.0*Y*Y-2.0*Z*Z) + +(2.0*X*Y+2.0*Z*W)*t1+1.0) / t2; + + s->b= ((2.0*X*Y+2.0*Z*W)*t1*(2.0*Y*Z+2.0*X*W) + -(1.0-2.0*Y*Y-2.0*Z*Z)*(2.0*Y*Z+2.0*X*W) + +2.0*Y*Z+2.0*X*W + -(2.0*X*Z-2.0*Y*W)*t1*(1.0-2.0*X*X-2.0*Z*Z) + +(2.0*X*Z-2.0*Y*W)*t1) + /(t1*t2); + + s->c= -(-(2.0*X*Z-2.0*Y*W)*(2.0*X*Y-2.0*Z*W) + +(1.0-2.0*Y*Y-2.0*Z*Z)*(2.0*Y*Z+2.0*X*W)-2.0*Y*Z-2.0*X*W) / t1; + + s->d= -(2.0*X*Z-2.0*Y*W)*(2.0*X*Y-2.0*Z*W) + +(1.0-2.0*Y*Y-2.0*Z*Z)*(2.0*Y*Z+2.0*X*W); + + s->e= 2.0*X*Y-2.0*Z*W; + + s->f= -((1.0-2.0*X*X-2.0*Y*Y)*t1-(2.0*Y*Z+2.0*X*W)*(2.0*X*Z+2.0*Y*W) + -2.0*X*Y+2.0*Z*W) / t2; + + s->g= ((1.0-2.0*X*X-2.0*Y*Y)*t1+t2*(2.0*X*Z+2.0*Y*W) + -(2.0*Y*Z+2.0*X*W)*(2.0*X*Z+2.0*Y*W)-t1) / (t1*t2); + + s->h= (-2.0*Y*Y-2.0*Z*Z)/t1; + if (debug) { + fprintf(stderr,"Result of decomp: a=%f, b=%f, c=%f, d= %f\n", + s->a,s->b,s->c,s->d); + fprintf(stderr," e= %f, f=%f, g=%f, h=%f\n", + s->e,s->f,s->g,s->h); + } + } + } + + static void trans_shear_adjust(ShearParams* s, TransParams* t_in, + TransParams* t_out) + { + t_out->dx= t_in->dx; + t_out->dy= t_in->dy - s->b*t_in->dx - s->a*t_in->dz; + t_out->dz= t_in->dz - s->c*t_in->dx; + } + + static int safe_quat( Quat* q ) { + /* Returns true if we are well away from singularities in the + * demoninators of the shear coefficient equations. + */ + double X,Y,Z,W,t1,t2,s1,s2; + + X= q->x; + Y= q->y; + Z= q->z; + W= q->w; + + /* Identity quaternion is safe */ + if (X==0 && Y==0 && Z==0) return 1; + else { + + t1= (2.0*X*Y-2.0*Z*W); + t2= (-(2.0*X*Z-2.0*Y*W)*t1 + +(1.0-2.0*Y*Y-2.0*Z*Z)*(2.0*Y*Z+2.0*X*W)); + + #define SQUARE(x) (x)*(x) + s1= sqrt(SQUARE(2.0*X*Y) + SQUARE(2.0*Z*W)); + s2= sqrt(SQUARE(2.0*X*Z*t1) + SQUARE(2.0*Y*W*t1) + + SQUARE(2.0*Y*Z) + SQUARE(2.0*X*W) + + SQUARE(4.0*Y*Y*Y*Z) + SQUARE(4.0*Y*Y*X*W) + + SQUARE(4.0*Y*Z*Z*Z) + SQUARE(4.0*X*Z*Z*W)); + #undef SQUARE + + if (debug) fprintf(stderr,"safe_quat: ratios %f, %f\n",t1/s1, t2/s2); + + return( (fabs(t1/s1) >= CANCELLATION_TOL) + && (fabs(t2/s2) >= CANCELLATION_TOL) ); + } + } + + static void step_repairing_quat( Quat* result, Quat* q, int step ) { + Quat qstep; + Quat tmp; + + if (debug) fprintf(stderr,"Repair step %d\n",step); + + switch (step % 3) { + case 0: quat_create(&qstep, REPAIR_QUAT_STEPSIZE, 0.0, 0.0, 1); break; + case 1: quat_create(&qstep, 0.0, REPAIR_QUAT_STEPSIZE, 0.0, 1); break; + case 2: quat_create(&qstep, 0.0, 0.0, REPAIR_QUAT_STEPSIZE, 1); break; + } + + quat_mult_right(result, &qstep); + if (debug) fprintf(stderr,"Trying (%f %f %f %f) with (%f %f %f %f)\n", + result->x, result->y, result->z, result->w, + q->x, q->y, q->z, q->w); + quat_copy(&tmp, result); + quat_mult_right(&tmp, q); + if (safe_quat(&tmp)) return; + else if (stepx,q->y,q->z,q->w); + } + + static void find_repairing_quat( Quat* result, Quat* q ) { + quat_identity(result); + step_repairing_quat( result, q, 0 ); + } + + static void shear_x( FComplex* image, double a, double b, double delta, + long nx, long ny, long nz, + double length_x, double length_y, double length_z ) + { + FComplex* here; + long halfx= nx/2; + long halfy= ny/2; + long halfz= nz/2; + int i; + int j; + int k; + double theta; + double x_phase_scale; + double y_phase_scale; + double z_phase_scale; + + /* Return if there is nothing to do */ + if (a==0.0 && b==0.0 && delta==0.0) return; + + if (debug) fprintf(stderr,"shear_x %f %f %f\n",a,b,delta); + + /* Step counter */ + count_shear_x++; + + /* FFT in x */ + fft3d( image, nx, ny, nz, +1, "x" ); + + /* delta is in voxels, a and b in fractional shears + * (typical range -1 to 1) + */ + x_phase_scale= 2.0*M_PI*delta/((double)nx); + y_phase_scale= 2.0*M_PI*a*length_y/(length_x*ny); + z_phase_scale= 2.0*M_PI*b*length_z/(length_x*nz); + + /* Apply phase changes */ + for (i=0; ireal; + t.imag= here->imag; + theta= (i-halfx)*x_phase_scale + + (j-halfy)*(i-halfx)*y_phase_scale + + (k-halfz)*(i-halfx)*z_phase_scale; + c= cos(theta); + s= sin(theta); + here->real= c*t.real - s*t.imag; + here->imag= c*t.imag + s*t.real; + } + } + } + + /* FFT back in x */ + fft3d( image, nx, ny, nz, -1, "x" ); + } + + static void shear_y( FComplex* image, double a, double b, double delta, + long nx, long ny, long nz, + double length_x, double length_y, double length_z ) + { + FComplex* here; + long halfx= nx/2; + long halfy= ny/2; + long halfz= nz/2; + int i; + int j; + int k; + double theta; + double x_phase_scale; + double y_phase_scale; + double z_phase_scale; + + /* Return if there is nothing to do */ + if (a==0.0 && b==0.0 && delta==0.0) return; + + if (debug) fprintf(stderr,"shear_y %f %f %f\n",a,b,delta); + + /* Step counter */ + count_shear_y++; + + /* FFT in y */ + fft3d( image, nx, ny, nz, +1, "y" ); + + /* delta is in voxels, a and b in fractional shears + * (typical range -1 to 1) + */ + y_phase_scale= 2.0*M_PI*delta/((double)ny); + z_phase_scale= 2.0*M_PI*a*length_z/(length_y*nz); + x_phase_scale= 2.0*M_PI*b*length_x/(length_y*nx); + + /* Apply phase changes */ + for (i=0; ireal; + t.imag= here->imag; + theta= (j-halfy)*y_phase_scale + + (k-halfz)*(j-halfy)*z_phase_scale + + (i-halfx)*(j-halfy)*x_phase_scale; + c= cos(theta); + s= sin(theta); + here->real= c*t.real - s*t.imag; + here->imag= c*t.imag + s*t.real; + } + } + } + + /* FFT back in y */ + fft3d( image, nx, ny, nz, -1, "y" ); + } + + static void shear_z( FComplex* image, double a, double b, double delta, + long nx, long ny, long nz, + double length_x, double length_y, double length_z ) + { + FComplex* here; + long halfx= nx/2; + long halfy= ny/2; + long halfz= nz/2; + int i; + int j; + int k; + double theta; + double x_phase_scale; + double y_phase_scale; + double z_phase_scale; + + /* Return if there is nothing to do */ + if (a==0.0 && b==0.0 && delta==0.0) return; + + if (debug) fprintf(stderr,"shear_z %f %f %f\n",a,b,delta); + + /* Step counter */ + count_shear_z++; + + /* FFT in z */ + fft3d( image, nx, ny, nz, +1, "z" ); + + /* delta is in voxels, a and b in fractional shears + * (typical range -1 to 1) + */ + z_phase_scale= 2.0*M_PI*delta/((double)nz); + x_phase_scale= 2.0*M_PI*a*length_x/(length_z*nx); + y_phase_scale= 2.0*M_PI*b*length_y/(length_z*ny); + + /* Apply phase changes */ + for (i=0; ireal; + t.imag= here->imag; + theta= (k-halfz)*z_phase_scale + + (i-halfx)*(k-halfz)*x_phase_scale + + (j-halfy)*(k-halfz)*y_phase_scale; + c= cos(theta); + s= sin(theta); + here->real= c*t.real - s*t.imag; + here->imag= c*t.imag + s*t.real; + } + } + } + + /* FFT back in z */ + fft3d( image, nx, ny, nz, -1, "z" ); + } + + void shift_only( TransParams* t, FComplex* moved_image, + long nx, long ny, long nz, + double length_x, double length_y, double length_z ) + { + /* Easy as pie, since pure shifts commute */ + shear_x(moved_image, 0.0, 0.0, t->dx, + nx, ny, nz, length_x, length_y, length_z); + shear_y(moved_image, 0.0, 0.0, t->dy, + nx, ny, nz, length_x, length_y, length_z); + shear_z(moved_image, 0.0, 0.0, t->dz, + nx, ny, nz, length_x, length_y, length_z); + + } + + void rot_x( Quat* q, TransParams* t, FComplex* moved_image, + long nx, long ny, long nz, + double length_x, double length_y, double length_z ) + { + TransParams adjusted_trans; + double alpha; + double delta; + + alpha= -(q->x/q->w); + delta= 2.0*q->x*q->w; + adjusted_trans.dx= t->dx; + adjusted_trans.dy= t->dy - alpha*t->dz; + adjusted_trans.dz= t->dz; + + shear_y(moved_image, alpha, 0.0, 0.0, nx, ny, nz, + length_x, length_y, length_z); + shear_z(moved_image, 0.0, delta, adjusted_trans.dz, nx, ny, nz, + length_x, length_y, length_z); + shear_y(moved_image, alpha, 0.0, adjusted_trans.dy, nx, ny, nz, + length_x, length_y, length_z); + shear_x(moved_image, 0.0, 0.0, adjusted_trans.dx, nx, ny, nz, + length_x, length_y, length_z); + } + + void rot_y( Quat* q, TransParams* t, FComplex* moved_image, + long nx, long ny, long nz, + double length_x, double length_y, double length_z ) + { + TransParams adjusted_trans; + double alpha; + double delta; + + alpha= -(q->y/q->w); + delta= 2.0*q->y*q->w; + adjusted_trans.dx= t->dx; + adjusted_trans.dy= t->dy; + adjusted_trans.dz= t->dz - alpha*t->dx; + + shear_z(moved_image, alpha, 0.0, 0.0, nx, ny, nz, + length_x, length_y, length_z); + shear_x(moved_image, 0.0, delta, adjusted_trans.dx, nx, ny, nz, + length_x, length_y, length_z); + shear_z(moved_image, alpha, 0.0, adjusted_trans.dz, nx, ny, nz, + length_x, length_y, length_z); + shear_y(moved_image, 0.0, 0.0, adjusted_trans.dy, nx, ny, nz, + length_x, length_y, length_z); + } + + void rot_z( Quat* q, TransParams* t, FComplex* moved_image, + long nx, long ny, long nz, + double length_x, double length_y, double length_z ) + { + TransParams adjusted_trans; + double alpha; + double delta; + + alpha= -(q->z/q->w); + delta= 2.0*q->z*q->w; + adjusted_trans.dx= t->dx - alpha*t->dy; + adjusted_trans.dy= t->dy; + adjusted_trans.dz= t->dz; + + shear_x(moved_image, alpha, 0.0, 0.0, nx, ny, nz, + length_x, length_y, length_z); + shear_y(moved_image, 0.0, delta, adjusted_trans.dy, nx, ny, nz, + length_x, length_y, length_z); + shear_x(moved_image, alpha, 0.0, adjusted_trans.dx, nx, ny, nz, + length_x, length_y, length_z); + shear_z(moved_image, 0.0, 0.0, adjusted_trans.dz, nx, ny, nz, + length_x, length_y, length_z); + } + + /* This routine applies appropriate phase shifts to the given + * image to cause it to be shifted by the given distances the + * next time it undergoes a full 3D FFT. If kspace_flag is + * set, the coming FFT is expected to be an inverse FFT; + * otherwise a forward FFT is assumed. + */ + void fshrot3d_set_shift_phases( double dx, double dy, double dz, + FComplex* image, + long nx, long ny, long nz, double length_x, + double length_y, double length_z, + int kspace_flag ) + { + FComplex* here; + long halfx= nx/2; + long halfy= ny/2; + long halfz= nz/2; + int i; + int j; + int k; + double theta; + double x_phase_scale; + double y_phase_scale; + double z_phase_scale; + + /* Return if there is nothing to do */ + if (dx==0.0 && dy==0.0 && dz==0.0) return; + + if (debug) fprintf(stderr,"fshrot3d_set_shift_phases %f %f %f\n", + dx, dy, dz); + + /* Step counter */ + count_set_phase++; + + /* No need for FFT here; we are assuming k-space data. */ + + /* dx, dy, dz are in voxels. Note that we need to reverse + * phases if we are doing the unexpected case (the future + * FFT is in the positive direction). + */ + x_phase_scale= 2.0*M_PI*dx/((double)nx); + y_phase_scale= 2.0*M_PI*dy/((double)ny); + z_phase_scale= 2.0*M_PI*dz/((double)nz); + if (!kspace_flag) { + x_phase_scale= -x_phase_scale; + y_phase_scale= -y_phase_scale; + /* No Z case because we are doing Z fft ourselves */ + } + + /* Apply phase changes */ + for (i=0; ireal; + t.imag= here->imag; + theta= (i-halfx)*x_phase_scale + + (j-halfy)*y_phase_scale + + (k-halfz)*z_phase_scale; + c= cos(theta); + s= sin(theta); + here->real= c*t.real - s*t.imag; + here->imag= c*t.imag + s*t.real; + } + } + } + + /* No need for FFT here; we are assuming k-space data. */ + + } + + /* In this implementation, we will use the decomposition: + * R = Sy(a,b)Sz(c,d)Sx(e,f)Sy(g,h) + * where R is the rotation matrix associated with q. Input and + * output data are assumed to be ordered such that x is fastest + * in memory. (Note that other expansions are used for special + * cases). + */ + void fourier_shift_rot3d( Quat* q, double dx, double dy, double dz, + FComplex* orig_image, + FComplex* moved_image, + long nx, long ny, long nz, + double length_x, double length_y, double length_z ) + { + TransParams trans; + int i; + + trans.dx= dx; + trans.dy= dy; + trans.dz= dz; + + /* Step counter */ + count_calls++; + + /* Copy into output buffer */ + for (i=0; ix==0.0 && q->y==0.0 && q->z==0.0) { + shift_only( &trans, moved_image, nx, ny, nz, + length_x, length_y, length_z ); + } + else if (q->x==0.0 && q->y==0.0 && fabs(q->w)>=CANCELLATION_TOL) { + rot_z( q, &trans, moved_image, nx, ny, nz, + length_x, length_y, length_z ); + } + else if (q->y==0.0 && q->z==0.0 && fabs(q->w)>=CANCELLATION_TOL) { + rot_x( q, &trans, moved_image, nx, ny, nz, + length_x, length_y, length_z ); + } + else if (q->z==0.0 && q->x==0.0 && fabs(q->w)>=CANCELLATION_TOL) { + rot_y( q, &trans, moved_image, nx, ny, nz, + length_x, length_y, length_z ); + } + else if (safe_quat(q)) { /* Check for numerical problems + * (singular points in shear coeffs). */ + + /* No problems; rotate and translate in 4 shears. */ + ShearParams q_params; + TransParams adjusted_trans; + + /* Find all the right parameters */ + quat_to_shears(q, &q_params); + trans_shear_adjust(&q_params, &trans, &adjusted_trans); + if (debug) fprintf(stderr,"Adjusted trans params %f %f %f\n", + adjusted_trans.dx, adjusted_trans.dy, + adjusted_trans.dz); + + /* Do the shears in place */ + shear_y(moved_image, q_params.g, q_params.h, 0.0, nx, ny, nz, + length_x, length_y, length_z); + shear_x(moved_image, q_params.e, q_params.f, adjusted_trans.dx, + nx, ny, nz, length_x, length_y, length_z); + shear_z(moved_image, q_params.c, q_params.d, adjusted_trans.dz, + nx, ny, nz, length_x, length_y, length_z); + shear_y(moved_image, q_params.a, q_params.b, adjusted_trans.dy, + nx, ny, nz, length_x, length_y, length_z); + } + else { + /* + * Found trouble, so generate a repairing quaternion and its inverse, + * and frame the shift and rotate with them. This results in 7 shears. + */ + Quat repairing_quat; + Quat repaired_q; + Quat repairing_conj; + ShearParams repaired_q_params; + ShearParams repairing_conj_params; + TransParams adjusted_trans; + + if (debug) fprintf(stderr,"Repair is needed\n"); + + /* Find the repairing quaternion, and mix it in with q */ + find_repairing_quat(&repairing_quat, q); + quat_copy(&repaired_q,&repairing_quat); + quat_mult_right(&repaired_q,q); + quat_copy(&repairing_conj, &repairing_quat); + quat_conjugate(&repairing_conj); + + /* Find all the right parameters */ + quat_to_shears(&repaired_q, &repaired_q_params); + quat_to_shears(&repairing_conj, &repairing_conj_params); + trans_shear_adjust(&repairing_conj_params, &trans, &adjusted_trans); + + /* Do the shears in place */ + shear_y(moved_image, repaired_q_params.g, repaired_q_params.h, 0.0, + nx, ny, nz, length_x, length_y, length_z); + shear_x(moved_image, repaired_q_params.e, repaired_q_params.f, 0.0, + nx, ny, nz, length_x, length_y, length_z); + shear_z(moved_image, repaired_q_params.c, repaired_q_params.d, 0.0, + nx, ny, nz, length_x, length_y, length_z); + shear_y(moved_image, + repairing_conj_params.g + repaired_q_params.a, + repairing_conj_params.h + repaired_q_params.b, + 0.0, nx, ny, nz, length_x, length_y, length_z); + shear_x(moved_image, + repairing_conj_params.e, repairing_conj_params.f, + adjusted_trans.dx, nx, ny, nz, length_x, length_y, length_z); + shear_z(moved_image, + repairing_conj_params.c, repairing_conj_params.d, + adjusted_trans.dz, nx, ny, nz, length_x, length_y, length_z); + shear_y(moved_image, + repairing_conj_params.a, repairing_conj_params.b, + adjusted_trans.dy, nx, ny, nz, length_x, length_y, length_z); + } + } + *** src/fmri/fshrot3d.h Wed Mar 24 17:24:04 1999 --- /SLOTH-U2/welling/fiasco/development/Fiasco4.2/src/fmri//fshrot3d.h Fri Feb 26 18:32:46 1999 *************** *** 0 **** --- 1,49 ---- + /************************************************************ + * * + * fshrot3d.h * + * * + * Permission is hereby granted to any individual or * + * institution for use, copying, or redistribution of * + * this code and associated documentation, provided * + * that such code and documentation are not sold for * + * profit and the following copyright notice is retained * + * in the code and documentation: * + * * + * Original programming by Joel Welling 2/99 * + ************************************************************/ + /* Header file for fshrot3d.c */ + + /* + Quat* q and the dx, dy, dz values specify rotation and shift (in voxels) + orig_image is input + moved_image is output + nx, ny, nz are grid dimensions + length_x, length_y, length_z are voxel edge lengths + */ + void fourier_shift_rot3d( Quat* q, double dx, double dy, double dz, + FComplex* orig_image, + FComplex* moved_image, + long nx, long ny, long nz, + double length_x, double length_y, double length_z); + + /* This routine simply sets phases in an image in such a way that + * a subsequent full 3D FFT causes a translation given by dx, dy, + * and dz (which are in voxels). The image is phase shifted in + * place. kspace_flag is nonzero if this is a k-space image, + * so that the expected coming FFT is in the inverse direction. + * / + void fshrot3d_set_shift_phases( double dx, double dy, double dz, + FComplex* image, + long nx, long ny, long nz, double length_x, + double length_y, double length_z, + int kspace_flag ); + + /* Turn debugging output on or off */ + void fshrot3d_set_debug( int flag ); + + /* Clear and get the counters for the shear operations (for diagnostics) */ + void fshrot3d_clear_shear_counts(); + void fshrot3d_get_shear_counts( int* xcount, int* ycount, int* zcount, + int* phase_count, + double* mean_shears_per_call ); + *** src/fmri/fshrot3d_help.help Wed Mar 24 17:24:04 1999 --- /SLOTH-U2/welling/fiasco/development/Fiasco4.2/src/fmri//fshrot3d_help.help Fri Feb 26 18:25:37 1999 *************** *** 0 **** --- 1,59 ---- + *FSHROT3DCalculation + + Fourier Shifts and rotations in 3D are implemented using the + "fshrot3d" package. This package decomposes a requested + rotation and translation into three or more shear operations + in 3D, and implements those shears using 1D Fourier transforms. + The elementary operations implement a shear, plus a translation + in the shear direction. + + If the operation is a simple translation (no rotation component), + three shearing operations implement it- one in each direction. + Hereafter rotations R may include a translation T. Pure rotations + will be denoted by lowercase r, such that + + R = T r + + + If an axis-aligned rotation is requested, three shears are used. + A fourth shear adds a translation along the rotational axis if + necessary. + Specifically, + + Rz = Sz Sy Sx Sy + + where Rz is a rotation about Z, Sx is a shear in the X direction, + and etc. Rotations about other axes can be found by permutation + of coordinates. + + A non-axis-aligned rotation can be implemented as four shears. + This package uses the expansion: + + R = Sy Sz Sx Sy + + The necessary coefficients for these shears are singular or + near-singular in some regions, however. If this situation is + found to occur, a correcting (pure) rotation rc is found which + shifts the rotation out of the bad region, and the rotation is + implemented as follows (bars denote the conjugate rotation): + + R = T r + __ + = T rc rc r + + = Sy Sz Sx Sy Sz Sx Sy + __ + by decomposing ( rc r ) and ( T rc ) each into four shears and + merging the first and last shears of each group. + + Note that the Fourier methods used here assume that the volume + being treated is periodic. As a result of this, if any shear + carries the "interesting" part of the data beyond the edge of + the volume, that data will reappear on the opposite edge and + produce strange artifacts. Thus this method is best used with + data confined to the region near the center of the volume, and + rotation and translation parameters that do not carry that + data to the boundary. + + m4include(../fmri/fft3d_help.help) + *** src/fmri/quat_tester.c Wed Mar 24 17:24:04 1999 --- /SLOTH-U2/welling/fiasco/development/Fiasco4.2/src/fmri//quat_tester.c Wed Feb 10 19:25:02 1999 *************** *** 0 **** --- 1,66 ---- + /************************************************************ + * * + * quat_tester.c * + * * + * Permission is hereby granted to any individual or * + * institution for use, copying, or redistribution of * + * this code and associated documentation, provided * + * that such code and documentation are not sold for * + * profit and the following copyright notice is retained * + * in the code and documentation: * + * * + * Original programming by Joel Welling 2/99 * + ************************************************************/ + /* Header file for quaternions and associated transforms */ + + #include + #include + #include "quaternion.h" + + #define DEG2RAD (M_PI/180.0) + + main( int argc, char* argv[] ) + { + Quat q; + Quat q_conj; + Transform t; + + if (argc != 1) { + fprintf(stderr,"Usage: %s\n",argv[0]); + fprintf(stderr," You will be prompted for Euler angles RzRyRx.\n"); + fprintf(stderr, + " An empty line causes previously given quats to be be multiplied.\n"); + exit(-1); + } + + quat_identity(&q); + + while (1) { + double tx, ty, tz; + int count; + char inbuf[81]; + printf("Enter 3 Euler angles X, Y, Z (used as RzRyRx) in degrees:\n"); + fgets(inbuf,80,stdin); + if (strlen(inbuf)<2) break; + else if ((count=sscanf(inbuf,"%lf %lf %lf\n",&tx,&ty,&tz)) == 3) { + Quat r; + quat_from_euler_RzRyRx(&r,DEG2RAD*tx,DEG2RAD*ty,DEG2RAD*tz); + quat_mult_right(&q,&r); + } + else printf("Try again.\n"); + } + + printf("Resulting quaternion: (%f %f %f %f)\n",q.x,q.y,q.z,q.w); + printf("Equivalent rotation matrix:\n"); + quat_to_trans(t,&q,0.0,0.0,0.0); + trans_dump(stdout,t); + trans_to_quat(&q_conj,t); + printf("Changing back to quaternion: (%f %f %f %f)\n",q_conj.x,q_conj.y, + q_conj.z,q_conj.w); + quat_conjugate(&q_conj); + quat_mult_right(&q,&q_conj); + printf("After multiplication by conjugate: get (%f %f %f %f)\n", + q.x,q.y,q.z,q.w); + quat_to_trans(t,&q,0.0,0.0,0.0); + trans_dump(stdout,t); + } *** src/fmri/quaternion.c Wed Mar 24 17:24:04 1999 --- /SLOTH-U2/welling/fiasco/development/Fiasco4.2/src/fmri//quaternion.c Wed Feb 17 20:07:32 1999 *************** *** 0 **** --- 1,242 ---- + /************************************************************ + * * + * quaternion.c * + * * + * Permission is hereby granted to any individual or * + * institution for use, copying, or redistribution of * + * this code and associated documentation, provided * + * that such code and documentation are not sold for * + * profit and the following copyright notice is retained * + * in the code and documentation: * + * * + * Original programming by Joel Welling 2/99 * + ************************************************************/ + /* Methods for quaternions and associated transforms */ + + #include + #include + #include "misc.h" + #include "quaternion.h" + + static char rcsid[] = "$Id: quaternion.c,v 1.3 1999/02/18 01:07:31 welling Exp $"; + + /* This package provides manipulations of unit quaternions. + Much code was taken from the "Matrix and Quaternions Faq". + Unfortunately the trans_to_quat algorithm given there seems + totally screwed up; this one was rederived using Maple. + */ + + static Quat* quat_alloc() + { + Quat* result; + if (!(result= (Quat*)malloc(sizeof(Quat)))) + Abort("Unable to allocate %d bytes for quaternion!\n",sizeof(Quat)); + return result; + } + + Quat* trans_to_quat( Quat* q_out, Transform t ) + { + double trace; + + if (!q_out) q_out= quat_alloc(); + + trace= t[0]+t[5]+t[10]+1.0; + if (trace>0.0) { + double s; + s= 0.5/sqrt(trace); + q_out->w= 0.25/s; + q_out->x = ( t[9] - t[6] ) * s; + q_out->y = ( t[2] - t[8] ) * s; + q_out->z = ( t[4] - t[1] ) * s; + } + else { + if (t[0]>t[5]) { + if (t[0]>t[10]) { + /* column 0 has greatest diagonal element */ + double S = 2.0/sqrt( 1.0 + t[0] - t[5] - t[10] ); /* = 1/X */ + q_out->x = 1.0 / S; + q_out->y = 0.25*(t[1] + t[4] )*S; + q_out->z = 0.25*(t[2] + t[8] )*S; + q_out->w = 0.25*(t[9] - t[6] )*S; + } + else { + /* column 2 has greatest diagonal element */ + double S = 2.0/sqrt( 1.0 + t[10] - t[0] - t[5] ); /* = 1/Z */ + q_out->x = 0.25*(t[2] + t[8] )*S; + q_out->y = 0.25*(t[6] + t[9] )*S; + q_out->z = 1.0 / S; + q_out->w = (t[4] - t[1] ) * S; + } + } + else { + if (t[5]>t[10]) { + /* column 1 has greatest diagonal element */ + double S = 2.0/sqrt( 1.0 + t[5] - t[0] - t[10] ); /* = 1/Y */ + q_out->x = 0.25*(t[1] + t[4])*S; + q_out->y = 1.0 / S; + q_out->z = 0.25*(t[6] + t[9])*S; + q_out->w = 0.25*(t[2] - t[8])*S; + } + else { + /* column 2 has greatest diagonal element */ + double S = 2.0/sqrt( 1.0 + t[10] - t[0] - t[5] ); /* = 1/Z */ + q_out->x = 0.25*(t[2] + t[8] )*S; + q_out->y = 0.25*(t[6] + t[9] )*S; + q_out->z = 1.0 / S; + q_out->w = (t[4] - t[1] ) * S; + } + } + } + + /* quat_normalize(q_out); */ + + return q_out; + } + + void quat_to_trans( Transform t_out, + Quat* q, double dx, double dy, double dz ) + { + double xx, xy, xz, xw, yy, yz, yw, zz, zw; + + xx = q->x * q->x; + xy = q->x * q->y; + xz = q->x * q->z; + xw = q->x * q->w; + + yy = q->y * q->y; + yz = q->y * q->z; + yw = q->y * q->w; + + zz = q->z * q->z; + zw = q->z * q->w; + + t_out[0] = 1 - 2 * ( yy + zz ); + t_out[1] = 2 * ( xy - zw ); + t_out[2] = 2 * ( xz + yw ); + + t_out[4] = 2 * ( xy + zw ); + t_out[5] = 1 - 2 * ( xx + zz ); + t_out[6] = 2 * ( yz - xw ); + + t_out[8] = 2 * ( xz - yw ); + t_out[9] = 2 * ( yz + xw ); + t_out[10] = 1 - 2 * ( xx + yy ); + + t_out[3] = dx; + t_out[7] = dy; + t_out[11]= dz; + t_out[12] = t_out[13] = t_out[14] = 0; + t_out[15] = 1; + } + + Quat* quat_copy( Quat* q_out, Quat* q_in ) + { + if (!q_out) q_out= quat_alloc(); + + q_out->x= q_in->x; + q_out->y= q_in->y; + q_out->z= q_in->z; + q_out->w= q_in->w; + + return q_out; + } + + double quat_magnitude( Quat* q ) + { + return( sqrt(q->w*q->w+q->x*q->x+ q->y*q->y+q->z*q->z) ); + } + + Quat* quat_normalize( Quat* q ) + { + double mag= quat_magnitude(q); + q->x /= mag; + q->y /= mag; + q->z /= mag; + q->w /= mag; + } + + Quat* quat_mult_right( Quat* q, Quat* f ) + { + Quat t; + t.w= q->w*f->w - (q->x*f->x + q->y*f->y + q->z*f->z); + t.x= q->y*f->z - q->z*f->y + q->w*f->x + q->x*f->w; + t.y= q->z*f->x - q->x*f->z + q->w*f->y + q->y*f->w; + t.z= q->x*f->y - q->y*f->x + q->w*f->z + q->z*f->w; + quat_copy(q,&t); + return q; + } + + Quat* quat_identity( Quat* q ) + { + if (!q) q= quat_alloc(); + q->w= 1.0; + q->x= q->y= q->z= 0.0; + return q; + } + + Quat* quat_conjugate( Quat* q ) + { + q->x= -q->x; + q->y= -q->y; + q->z= -q->z; + return q; + } + + Quat* quat_create( Quat* q, double x, double y, double z, int w_pos ) + { + if (!q) q= quat_alloc(); + q->x= x; + q->y= y; + q->z= z; + q->w= sqrt(1.0 - (x*x + y*y + z*z)); + if (!w_pos) q->w= -q->w; + return q; + } + + Quat* quat_from_axis_angle( Quat* q, double x, double y, double z, + double theta ) + { + double sin_a; + double cos_a; + + if (!q) q= quat_alloc(); + + sin_a= sin(0.5*theta); + cos_a= cos(0.5*theta); + + q->x= x * sin_a; + q->y= y * sin_a; + q->z= z * sin_a; + q->w= cos_a; + quat_normalize( q ); + + return q; + } + + Quat* quat_from_euler_RzRyRx( Quat* q, double x_angle, double y_angle, + double z_angle ) + { + Quat qx; + Quat qy; + Quat qz; + + if (!q) q= quat_alloc(); + + quat_from_axis_angle( &qx, 1.0, 0.0, 0.0, x_angle ); + quat_from_axis_angle( &qy, 0.0, 1.0, 0.0, y_angle ); + quat_from_axis_angle( &qz, 0.0, 0.0, 1.0, z_angle ); + + quat_copy(q, &qz); + quat_mult_right(q, &qy); + quat_mult_right(q, &qx); + + return q; + } + + void trans_dump( FILE* ofile, Transform t ) + { + fprintf(ofile," ( %5f %5f %5f %5f )\n",t[0],t[1],t[2],t[3]); + fprintf(ofile," ( %5f %5f %5f %5f )\n",t[4],t[5],t[6],t[7]); + fprintf(ofile," ( %5f %5f %5f %5f )\n",t[8],t[9],t[10],t[11]); + fprintf(ofile," ( %5f %5f %5f %5f )\n",t[12],t[13],t[14],t[15]); + } *** src/fmri/quaternion.h Wed Mar 24 17:24:04 1999 --- /SLOTH-U2/welling/fiasco/development/Fiasco4.2/src/fmri//quaternion.h Wed Feb 10 19:22:58 1999 *************** *** 0 **** --- 1,49 ---- + /************************************************************ + * * + * quaternion.h * + * * + * Permission is hereby granted to any individual or * + * institution for use, copying, or redistribution of * + * this code and associated documentation, provided * + * that such code and documentation are not sold for * + * profit and the following copyright notice is retained * + * in the code and documentation: * + * * + * Original programming by Joel Welling 2/99 * + ************************************************************/ + /* Header file for quaternions and associated transforms */ + + typedef struct quaternion_struct { + double x; + double y; + double z; + double w; + } Quat; /* a quaternion */ + + typedef double Transform[16]; /* a homogeneous transform */ + + void trans_dump( FILE* ofile, Transform t ); + + Quat* trans_to_quat( Quat* q_out, Transform t ); /* translations ignored */ + + void quat_to_trans( Transform t_out, + Quat* q, double dx, double dy, double dz ); + + Quat* quat_copy( Quat* q_out, Quat* q_in ); + + Quat* quat_normalize( Quat* q ); + + Quat* quat_mult_right( Quat* q, Quat* factor ); + + Quat* quat_identity( Quat* q ); + + Quat* quat_conjugate( Quat* q ); + + Quat* quat_create( Quat* q, double x, double y, double z, int w_pos ); + + Quat* quat_from_axis_angle( Quat* q, double x, double y, double z, + double theta ); + + Quat* quat_from_euler_RzRyRx( Quat* q, double x_angle, double y_angle, + double z_angle ); + *** src/mri_util/mri_pad.c Wed Mar 24 17:24:04 1999 --- /SLOTH-U2/welling/fiasco/development/Fiasco4.2/src/mri_util//mri_pad.c Thu Feb 11 19:56:50 1999 *************** *** 0 **** --- 1,362 ---- + /************************************************************ + * * + * mri_pad.c * + * * + * Permission is hereby granted to any individual or * + * institution for use, copying, or redistribution of * + * this code and associated documentation, provided * + * that such code and documentation are not sold for * + * profit and the following copyright notice is retained * + * in the code and documentation: * + * Copyright (c) 1997 Department of Statistics, * + * Carnegie Mellon University * + * * + * Original programming by Joel Welling, 5/98 * + ************************************************************/ + /************************************************************* + + DESCRIPTION OF MRI_PAD + + mri_pad takes a pgh MRI dataset of any type, and outputs + a dataset of the same type containing a pad of the data. + The pad represents a subrange of one of the dimensions of + the data. The pad operation acts on all chunks in the + dataset. + + mri_pad [-d dim] -r range [-o offset] [-f fillval] [-v] infile outfile + + -d dim (where dim is a single character) specifies the dimension + to be padded (default t). Valid values are + one character long; case is ignored. + -r range specifies the new extent of that dimension. This + value must be greater than or equal to the original extent + plus the offset. + -f number to use as padding (default 0). This value is read as + a float, then converted to match the needed data type(s). + -o offset specifies the value along the selected dimension at + which the original data begins. (default 0) + -v sets verbose output + infile specifies the input dataset + outfile specifies the output dataset + + **************************************************************/ + + #include + #include + #include + #include + #include "mri.h" + #include "fmri.h" + #include "misc.h" + #include "stdcrg.h" + + #define KEYBUF_SIZE 512 + + static char rcsid[] = "$Id: mri_pad.c,v 1.1 1999/02/12 00:56:30 welling Exp $"; + + static MRI_Dataset *Input = NULL, *Output = NULL; + static char selected_dim[512]= "t"; + static int offset; + static int range; + static int data_changed= 0; + static float fillval= 0.0; + static char* progname; + static int verbose_flg= 0; + + static void safe_copy(char* str1, char* str2) { + strncpy(str1, str2, KEYBUF_SIZE); + str1[KEYBUF_SIZE-1]= '\0'; + } + + static void safe_concat(char* str1, char* str2) { + strncat(str1, str2, (KEYBUF_SIZE-strlen(str1))-1); + } + + static int get_chunk_type(MRI_Dataset* ds, char* chunk) + { + char key_buf[KEYBUF_SIZE]; + + safe_copy(key_buf,chunk); + safe_concat(key_buf,".datatype"); + if (mri_has(ds, key_buf)) { + char* type_name= mri_get_string(ds, key_buf); + if (!strcmp(type_name,"uint8")) return MRI_UNSIGNED_CHAR; + else if (!strcmp(type_name,"int16")) return MRI_SHORT; + else if (!strcmp(type_name,"int32")) return MRI_INT; + else if (!strcmp(type_name,"float32")) return MRI_FLOAT; + else if (!strcmp(type_name,"float64")) return MRI_DOUBLE; + else Abort("%s: unknown data type for key %s!\n",progname,key_buf); + } + else Abort("%s: missing tag %s!\n",progname,key_buf); + + return 0; /* not reached */ + } + + static int safe_get_extent(MRI_Dataset* ds, char* chunk, char* dim) + { + char key_buf[KEYBUF_SIZE]; + char dim_buf[4]; + dim_buf[0]= *dim; + dim_buf[1]= '\0'; + safe_copy(key_buf,chunk); + safe_concat(key_buf,".extent."); + safe_concat(key_buf,dim_buf); + if (mri_has(Input,key_buf)) return mri_get_int(Input,key_buf); + else Abort("%s: input missing tag %s!\n",progname,key_buf); + return 0; /* not reached */ + } + + static void calc_sizes(char* this_chunk, char* dimstr, + int* fast_blocksize_out, int* slow_blocksize_out ) { + /* This routine will fail if selected_dim is not in dimstr! */ + int fast_blocksize; + int slow_blocksize; + char* this_dim; + + fast_blocksize= 1; + + this_dim= dimstr; + while (*this_dim != *selected_dim) + fast_blocksize *= safe_get_extent(Input,this_chunk,this_dim++); + + this_dim++; /* step over selected dim */ + + slow_blocksize= 1; + while (*this_dim) + slow_blocksize *= safe_get_extent(Input, this_chunk, this_dim++); + + *fast_blocksize_out= fast_blocksize; + *slow_blocksize_out= slow_blocksize; + } + + static void transfer_data(char* this_chunk, + int fast_blksize, int slow_blksize, + int selected_extent) + { + int in_offset= 0; + int out_offset= 0; + int ifast; + int islow; + int collective_blksize; + int out_blksize; + int type; + int typesize; + void* obuf; + void* ibuf; + int i; + + type= get_chunk_type(Input,this_chunk); + typesize= get_typesize(Input,this_chunk); + + collective_blksize= fast_blksize * selected_extent; + out_blksize= fast_blksize * range; + if (!(obuf= malloc(out_blksize*typesize))) { + Abort("%s: unable to allocate %d bytes!\n",progname, + out_blksize*typesize); + } + + for (islow=0; islow %d at %d\n", + type, collective_blksize, in_offset, + out_blksize, out_offset); + in_offset += collective_blksize; + out_offset += out_blksize; + } + + free(obuf); + } + + static void pad_chunk(char* this_chunk) { + char key_buf[KEYBUF_SIZE]; + char* dimstr; + int selected_extent; + + if (verbose_flg) fprintf(stderr,"padding chunk <%s>\n",this_chunk); + safe_copy(key_buf, this_chunk); + safe_concat(key_buf, ".dimensions"); + if (mri_has(Input,key_buf)) { + dimstr= mri_get_string(Input,key_buf); + if (verbose_flg) fprintf(stderr,"dimstr %s\n",dimstr); + if (strchr(dimstr,*selected_dim)) { + safe_copy(key_buf, this_chunk); + safe_concat(key_buf, ".extent."); + safe_concat(key_buf, selected_dim); + if (mri_has(Input,key_buf)) + selected_extent= mri_get_int(Input,key_buf); + else Abort("%s: input missing tag %s!\n",progname,key_buf); + if (selected_extent+offset>range) + Abort("%s: chunk <%s> is too wide to pad!\n",progname,this_chunk); + if (verbose_flg) + fprintf(stderr,"selected dim extent on input %d\n",selected_extent); + if (selected_extent 1 ) && !strcmp( argv[1], "-help" ) ) + { + if( argc == 2 ) + Help( "selecttopic" ); + else + Help( argv[2] ); + } + + /*** Parse command line ***/ + + cl_scan( argc, argv ); + + /* Get filenames */ + cl_get("d", "%option %s[t]",selected_dim); + if (strlen(selected_dim)>1) { + fprintf(stderr,"%s: Selected dim name must be 1 char long.\n",argv[0]); + Help("usage"); + exit(-1); + } + *selected_dim= tolower(*selected_dim); + if (!cl_get("r", "%option %d",&range)) { + fprintf(stderr,"%s: Index range not given.\n",argv[0]); + Help("usage"); + exit(-1); + } + cl_get("o", "%option %d[0]",&offset); + cl_get("f", "%option %f[0.0]",&fillval); + if (!cl_get("", "%s", infile)) { + fprintf(stderr,"%s: Input file name not given.\n",argv[0]); + Help( "usage" ); + exit(-1); + } + if (!cl_get("", "%s", outfile)) { + fprintf(stderr,"%s: Output file name not given.\n",argv[0]); + Help( "usage" ); + exit(-1); + } + verbose_flg= cl_present("v"); + if (cl_cleanup_check()) { + int i; + fprintf(stderr,"%s: invalid argument in command line:\n ",argv[0]); + for (i=0; i=KEYBUF_SIZE) + Abort("%s: key too long!\n",argv[0]); + if (!strchr(this_key,'.')) { + strncpy(this_chunk, this_key, KEYBUF_SIZE); + this_chunk[KEYBUF_SIZE-1]= '\0'; + pad_chunk(this_chunk); + } + } + + /* Write and close data-sets */ + mri_close_dataset( Input ); + mri_close_dataset( Output ); + + Message( "# Padding complete.\n" ); + if (!data_changed) + Message("# Warning: input and output datasets identical!\n"); + + return 0; + } + *** src/mri_util/mri_pad_help.help Wed Mar 24 17:24:04 1999 --- /SLOTH-U2/welling/fiasco/development/Fiasco4.2/src/mri_util//mri_pad_help.help Thu Feb 11 19:56:50 1999 *************** *** 0 **** --- 1,28 ---- + *Usage + + mri_pad takes a pgh MRI dataset of any type, and outputs a + dataset of the same type containing a padded version of the + data. This version contains the original data, with one + dimension extended by padding. The pad operation acts on + all chunks in the dataset. + + mri_pad [-d dim] -r range [-o offset] [-f fillval] [-v] infile outfile + + -d dim (where dim is a single character) specifies the dimension + to be padded (default t). Valid values are + one character long; case is ignored. + -r range specifies the new extent of that dimension. This + value must be greater than or equal to the original extent + plus the offset. + -f number to use as padding (default 0). This value is read as + a float, then converted to match the needed data type(s). + -o offset specifies the value along the selected dimension at + which the original data begins. (default 0) + -v sets verbose output + infile specifies the input dataset + outfile specifies the output dataset + + Chunks not containing the selected dimension are copied verbatim. + Any chunk containing too many entries to be padded as specified + causes mri_pad to abort. + *** src/ireg3d/Makefile Wed Mar 24 17:24:48 1999 --- /SLOTH-U2/welling/fiasco/development/Fiasco4.2/src/ireg3d//Makefile Wed Feb 17 20:03:46 1999 *************** *** 0 **** --- 1,36 ---- + # + # Makefile for ireg3d + # + # Copyright (c) 1996 Department of Statistics, Carnegie Mellon University + # + # HISTORY + # 12/96 Written by William F. Eddy (CMU) + # + + PKG = ireg3d + PKG_MAKEBINS = $(CB)/ireg3d + + PKG_LIBS = -lfmri -lmri -lpar -lbio -lacct -lmisc -lcrg -lm + + MAKEFILES= Makefile + CSOURCE= ireg3d.c + DOCFILES= ireg3d_help.help + + include ../Makefile.pkg + + LIBFILES= $L/libmri.a $L/libpar.a $L/libbio.a $L/libarray.a $L/libmisc.a \ + $L/libacct.a + HDRS= fmri.h mri.h par.h bio.h misc.h acct.h stdcrg.h + + $O/ireg3d.o: ireg3d.c + $(CC_RULE) + + $O/ireg3d_help.o: ireg3d_help.help + $(HELP_RULE) + + $(CB)/ireg3d: $O/ireg3d.o $O/ireg3d_help.o $(LIBFILES) + $(SINGLE_HELP_LD) + + releaseprep: + echo "no release prep from " `pwd` + *** src/ireg3d/depend.mk.dummy Wed Mar 24 17:24:48 1999 --- /SLOTH-U2/welling/fiasco/development/Fiasco4.2/src/ireg3d//depend.mk.dummy Wed Feb 10 19:27:03 1999 *************** *** 0 **** --- 1,2 ---- + # This is a harmless dummy dependencies file. + *** src/ireg3d/ireg3d.c Wed Mar 24 17:24:48 1999 --- /SLOTH-U2/welling/fiasco/development/Fiasco4.2/src/ireg3d//ireg3d.c Fri Feb 26 18:39:14 1999 *************** *** 0 **** --- 1,374 ---- + /************************************************************ + * * + * ireg.c * + * * + * Permission is hereby granted to any individual or * + * institution for use, copying, or redistribution of * + * this code and associated documentation, provided * + * that such code and documentation are not sold for * + * profit and the following copyright notice is retained * + * in the code and documentation: * + * Copyright (c) 1995 Department of Statistics, * + * Carnegie Mellon University * + * * + * Original programming by Mark Fitzgerald 2-95 * + * 5-96: Pittsburgh Format, Mark Fitzgerald * + ************************************************************/ + + #include + #include + #include + #include "mri.h" + #include "fmri.h" + #include "stdcrg.h" + #include "misc.h" + + static char rcsid[] = "$Id: ireg3d.c,v 1.4 1999/02/26 23:39:14 welling Exp $"; + + /* Access for 3D arrays */ + #define MEM_BCK(matrix,nx,ny,nz,x,y,z) matrix[((((z)*ny)+(y))*nx)+(x)] + #define MEM(matrix,nx,ny,nz,x,y,z) matrix[((((x)*ny)+(y))*nz)+(z)] + + typedef struct regpar3d_struct { + Quat q; + double x; + double y; + double z; + } RegPar3D; + + static int debug_flag= 0; + + static void load_reg_params(char* parfile, RegPar3D* par, long dt) + { + FILE *fp = NULL; + char scanline[512]; + long linenum, numread, tt; + Quat q; + double x_shift, y_shift, z_shift; + + fp = efopen( parfile, "r" ); + linenum = -1; + while( !feof( fp ) && !ferror( fp ) ) + { + linenum++; + + /* Scan a line, ignoring comments (which begin with '#') */ + if (fgets(scanline, sizeof(scanline), fp) + && strlen(scanline)>0 && scanline[0] != '#') { + numread = sscanf( scanline, "%ld%lf%lf%lf%lf%lf%lf%lf%*f", + &tt, &(q.x), &(q.y), &(q.z), &(q.w), + &x_shift, &y_shift, &z_shift ); + if( numread < 8 ) + { + Warning( 1, "Line %ld of %s is too short (%ld) -- Ignoring.\n", + linenum, parfile, numread ); + continue; + } + + /* Check to see if image and slice numbers are in bounds */ + if( ( tt < 0 ) || ( tt >= dt ) ) + { + Warning( 1, + "Image out of bounds (%ld) for line %ld of %s -- Ignoring.\n", + tt, linenum, parfile ); + continue; + } + quat_normalize(&q); + + /* Put parameters into appropriate storage */ + quat_copy(&(par[tt].q),&q); + par[tt].x= x_shift; + par[tt].y= y_shift; + par[tt].z= z_shift; + if (debug_flag) + fprintf(stderr,"%d: loaded (%f %f %f %f) %f %f %f\n", + tt,par[tt].q.x,par[tt].q.y,par[tt].q.z,par[tt].q.w, + par[tt].x,par[tt].y,par[tt].z); + } + } + efclose( fp ); + } + + main( argc, argv ) + int argc; + char **argv; + { + + MRI_Dataset *Input = NULL, *Output = NULL; + char infile[512], hdrfile[512], outfile[512], parfile[512]; + long dv, dx, dy, dz, dt; + FComplex *c_image = NULL, *c_corr_image = NULL, *c_image_in= NULL; + float *image_in = NULL, *corr_image = NULL; + long x, y, t, z; + RegPar3D* par= NULL; + Quat q; + long block_offset; + int block_size; + double length_x, length_y, length_z; + int shear_count_x, shear_count_y, shear_count_z, phase_count; + double shears_per_call; + double xvoxel, yvoxel, zvoxel; + + /* Print version number */ + Message( "# %s\n", rcsid ); + + /* Check to see if help was requested */ + if( ( argc > 1 ) && !strcmp( argv[1], "-help" ) ) + { + if( argc == 2 ) + Help( "selecttopic" ); + else + Help( argv[2] ); + } + + /*** Parse command line ***/ + + cl_scan( argc, argv ); + + /* Get filenames */ + cl_get( "dataout|d", "%option %s[%]", ".dat", outfile ); + cl_get( "headerout|h", "%option %s[%]", "ireg3d.mri", hdrfile ); + cl_get( "input|i", "%option %s[%]", "input.mri", infile ); + cl_get( "parameters|p", "%option %s[%]", "reg3d.par", parfile ); + if (!cl_get( "x|xvoxel", "%option %lf", &xvoxel )) { + fprintf(stderr,"%s: required argument xvoxel omitted.\n",argv[0]); + Help("usage"); + exit(-1); + } + if (!cl_get( "y|yvoxel", "%option %lf", &yvoxel )) { + fprintf(stderr,"%s: required argument yvoxel omitted.\n",argv[0]); + Help("usage"); + exit(-1); + } + if (!cl_get( "z|zvoxel", "%option %lf", &zvoxel )) { + fprintf(stderr,"%s: required argument zvoxel omitted.\n",argv[0]); + Help("usage"); + exit(-1); + } + debug_flag= cl_present("debug"); + + if (cl_cleanup_check()) { + int i; + fprintf(stderr,"%s: invalid argument in command line:\n ",argv[0]); + for (i=0; i 2 ) ) + Error( "%s takes only reals or complex numbers of the form (v)xyzt.", + argv[0] ); + } + else + Error( "%s does not have the images.dimensions key.", infile ); + + /* Set output dataset */ + Output = mri_copy_dataset( hdrfile, Input ); + mri_set_string( Output, "images.datatype", "float32" ); + mri_set_string( Output, "images.file", outfile ); + mri_set_string( Output, "images.dimensions", "vxyzt" ); + mri_set_int( Output, "images.extent.v", (int) dv ); + + /* Set parameters in local variables */ + if( !mri_has( Input, "images.extent.t" ) || + !mri_has( Input, "images.extent.x" ) || + !mri_has( Input, "images.extent.y" ) || + !mri_has( Input, "images.extent.z" ) ) + Error( "images.extent key(s) missing from header." ); + dt = mri_get_int( Input, "images.extent.t" ); + dx = mri_get_int( Input, "images.extent.x" ); + dy = mri_get_int( Input, "images.extent.y" ); + dz = mri_get_int( Input, "images.extent.z" ); + if( ( dt <= 0 ) || ( dx <= 0 ) || ( dy <= 0 ) || ( dz <= 0 ) ) + Error( "images.extent key(s) is non-positive." ); + + /* Allocate image and parameter storage */ + if (!(c_image= (FComplex*)malloc(dx*dy*dz*sizeof(FComplex)))) + Abort("%s: unable to allocate %d bytes!\n", + argv[0],dx*dy*dz*sizeof(FComplex)); + if (!(c_corr_image= (FComplex*)malloc(dx*dy*dz*sizeof(FComplex)))) + Abort("%s: unable to allocate %d bytes!\n", + argv[0],dx*dy*dz*sizeof(FComplex)); + if( dv == 1 ) { + if (!(corr_image= (float*)malloc(dx*dy*dz*sizeof(float)))) + Abort("%s: unable to allocate %d bytes!\n", + argv[0],dx*dy*dz*sizeof(float)); + } + if (!(par= (RegPar3D*)malloc(dt*sizeof(RegPar3D)))) + Abort("%s: unable to allocate %d bytes!\n", + dt*sizeof(Transform)); + + /* Initialize registration parameters */ + /* to zero in case they are missing */ + quat_identity(&q); + for( t = 0; t < dt; t++ ) { + quat_copy(&(par[t].q), &q); + par[t].x= par[t].y= par[t].z= 0.0; + } + + /* Read in estimated registration parameters */ + load_reg_params( parfile, par, dt ); + + /* If real-valued input, set imaginary part */ + /* of complex image to zero */ + if( dv == 1 ) { + for(x=0; x + #include + #include "mri.h" + #include "fmri.h" + #include "misc.h" + #include "stdcrg.h" + #include "polyfit.h" + + /* Empirically derived constants */ + #define DEFAULT_WINDOW 20 + #define DEFAULT_THRESHOLD 25.0 + + /* Global parameters */ + int window; + float threshold; + + int find_next_peak(float *, int, int); + void process_data(float *, float *, int); + + /* process_data: */ + /* does all the actual data processing taking the cardio data in the */ + /* input buffer and making our cooresponding ramping function in the */ + /* output buffer where n is the size of the buffers. */ + + void process_data(float *in_buffer, float *out_buffer, int n) { + float coeffs[3], slope; + int i, j, k, w2 = window/2; + + float *temp_buffer = (float *)malloc(n * sizeof(float)); + + /* Set the temp and output buffer contents to 0 */ + for(i = 0; i < n; i++) { + out_buffer[i] = 0; + temp_buffer[i] = 0; + } + + /* Make a first pass over the data fitting parabolas */ + for(i = w2; i < n - w2; i++) { + polyfitlin(in_buffer + i - w2, window, coeffs, 3); + temp_buffer[i] = coeffs[2] * 100; + } + + /* Now go through and find the peaks from the fitting stage */ + /* at each peak, we seek a local minimum in the orginal */ + /* data around it. */ + + for(i = w2; i < n - w2; i++) + if(temp_buffer[i-1] < temp_buffer[i] && + temp_buffer[i+1] < temp_buffer[i] && + temp_buffer[i] > threshold) { + j = i; + while((in_buffer[j] > in_buffer[j+1] || + in_buffer[j] > in_buffer[j-1]) && + (j > 1 && j < n-1)) + if(in_buffer[j+1] < in_buffer[j]) j = j+1; + else j = j-1; + out_buffer[j] = 1; + } + + /* Now linearize everything */ + + i = -1; + while((i = find_next_peak(out_buffer, i, n)) != -1) { + j = find_next_peak(out_buffer, i, n); + if(j == -1) break; + slope = 1.0/(j - i - 1); + for(k = i + 1; k < j; k++) + out_buffer[k] = (k - i - 1) * slope; + } + + /* Now we have to take care of the sections at the beginning and + end of the data that just didn't fit into our window for the + we copy the slope from the first and last sections where we were + able to identify peaks, and fill in the gaps. */ + + i = find_next_peak(out_buffer, -1, n); + j = find_next_peak(out_buffer, i, n); + for(k = 0; k < i; k++) out_buffer[k] = out_buffer[j - i + k]; + + for(i = n - 1; i >= 0 && out_buffer[i] == 0; i--); + for(j = i - 1; j >= 0 && out_buffer[j] != 1; j--); + for(k = i + 1; k < n; k++) out_buffer[k] = out_buffer[j - i + k]; + + } + + /* find_next_peak: */ + /* given the location of the previous peak, finds the next */ + /* peak. if none is to be found, returns -1 */ + + int find_next_peak(float *buffer, int start, int size) { + int i = start + 1; + while(i < size && buffer[i] != 1) i++; + if(i == size) return -1; + return i; + } + + int main(int argc, char *argv[]) { + MRI_Dataset *input, *output; + int chunk_size, num_elements, i, j, k; + float *input_buffer, *output_buffer; + char chunk[512]; + char infile[512]; + char outfile[512]; + char tag[600]; + + /* Check to see if help was requested */ + if( ( argc > 1 ) && !strcmp( argv[1], "-help" ) ) + { + if( argc == 2 ) + Help( "selecttopic" ); + else + Help( argv[2] ); + } + + /* First, parse the command line */ + + cl_scan(argc, argv); + + /* Then, get the options, input and output filenames */ + + cl_get("c", "%option %s[%]", "physio", chunk); + cl_get("w", "%option %d[%]", DEFAULT_WINDOW, &window); + cl_get("t", "%option %f[%]", DEFAULT_THRESHOLD, &threshold); + + if(!cl_get("", "%s", infile)) { + fprintf(stderr, "%s: Input file name not given.\n", argv[0]); + exit(-1); + } + if(!cl_get("", "%s", outfile)) { + fprintf(stderr, "%s: Output file name not given.\n", argv[0]); + exit(-1); + } + + if(cl_cleanup_check()) { + int i; + fprintf(stderr, "%s: invalid argument in command line:\n", argv[0]); + for(i = 0; i < argc; i++) fprintf(stderr, "%s ", argv[i]); + fprintf(stderr, "\n"); + Help("usage"); + exit(-1); + } + + /* Now that we are done the command-line parsing, lets copy the + whole dataset and open everything up so we can start messing + around with it. If the named chunk isn't in the dataset, + then we abort. We also abort if the input and output filenames + are the same */ + + if(!strcmp(infile, outfile)) { + fprintf(stderr, "%s: input and output files must be distinct.\n", argv[0]); + exit(-1); + } + input = mri_open_dataset(infile, MRI_READ); + if(!mri_has(input, chunk)) { + fprintf(stderr, "%s: requested chunk %s not present.\n", argv[0], chunk); + exit(-1); + } + strcpy(tag, chunk); + strcat(tag, ".dimensions"); + if (strcmp(mri_get_string(input,tag), "t")) { + fprintf(stderr,"%s: input chunk does not have dimensions 't'\n", argv[0]); + exit(-1); + } + output = mri_copy_dataset(outfile, input); + + /* Get the length of the chunk */ + strcpy(tag, chunk); + strcat(tag, ".extent.t"); + chunk_size = mri_get_int(input, tag); + + /* Now we read everything in and get the output buffer. Since the physio + data sets should not be ridiculously huge, everything gets read in + at the same time. */ + + input_buffer = mri_get_chunk(input, chunk, chunk_size, 0, MRI_FLOAT); + output_buffer = mri_get_chunk(output, chunk, chunk_size, 0, MRI_FLOAT); + + /* Then we process the data, write out the results and close our + files. Then there is much happiness and joy. */ + + process_data(input_buffer, output_buffer, chunk_size); + mri_set_chunk(output, chunk, chunk_size, 0, MRI_FLOAT, output_buffer); + + mri_close_dataset(input); + mri_close_dataset(output); + + exit(0); + } *** src/physio/cardiopeaks_help.help Wed Mar 24 17:25:21 1999 --- /SLOTH-U2/welling/fiasco/development/Fiasco4.2/src/physio//cardiopeaks_help.help Tue Feb 9 14:44:26 1999 *************** *** 0 **** --- 1,27 ---- + *Usage + + cardiopeaks takes a pgh MRI dataset containing a cardio chunk + and creates a new MRI file containing a chunk with a ramping + function cooresponding to the heart beats in the original file. + Specifically, the output file values run linearly from 0.0 + immediately following one beat up to 1.0 the moment before + the next. Estimates are made for the values before and after + the first and last peaks. The input chunk must be a pure time + series, that is, it's dimensions must be "t". + + Parameters: + + cardiopeaks [-c chunk] [-w window_size] [-t threshold] infile outfile + + -c chunk specifices the name of the chunk to use + (default "physio") + -w window_size specifies the size of the window to + use for filtering the data when looking for peaks. + (default 20 samples) + -t threshold specifies the threshold to use for + deciding whether a peak is significant or not + (default 25) + infile specifies the input dataset + outfile specifies the output dataset + + *** src/physio/polyfit.c Wed Mar 24 17:25:21 1999 --- /SLOTH-U2/welling/fiasco/development/Fiasco4.2/src/physio//polyfit.c Tue Feb 9 12:21:09 1999 *************** *** 0 **** --- 1,154 ---- + #include + #include + #include + + float polyfit(float *, int *, int, float *, int); + float polyfitlin(float *, int, float *, int); + void solve(double *, float *, int); + + /* polyfitlin: */ + /* Arguments - */ + /* data : an array of floats containing data points */ + /* n : # of data points */ + /* beta : array of floats for holding the coefficients */ + /* order : order of the polynomial to fit to */ + /* Returns - */ + /* mean squared error of the data with the fitted polynomial */ + + float polyfitlin(float *data, int n, float *beta, int order) { + int *x, i; + + x = (int *) malloc(n * sizeof(int)); + for(i = 0; i < n; i++) + x[i] = i; + + return polyfit(data, x, n, beta, order); + + } + + /* polyfit: */ + /* Arguments - */ + /* data : an array of float containing data points */ + /* x : an array of x values for the data points */ + /* n : # of data points */ + /* beta : array of floats for holding the cofficients */ + /* order : order of the polynomial to fit to */ + /* Returns - */ + /* mean squared error of the data with the fitted polynomial */ + + float polyfit(float *data, int *x, int n, float *beta, int order) { + int i, j, k, k2, k3; + double sum, xi, val, error; + + double *basis = (double *) malloc(order * n * sizeof(double)); + double *alpha = (double *) malloc(order * order * sizeof(double)); + + /* Calculate basis */ + + for(i = 0; i < order; i++) + for(j = 0; j < n; j++) { + k = i + j * order; + if(i == 0) basis[k] = 1; + else basis[k] = x[j] * basis[k - 1]; + } + + /* Calculate alpha */ + + for(i = 0; i < order; i++) { + for(j = 0; j <= i; j++) { + sum = 0; + for(k = 0; k < n; k++) { + k2 = i + k * order; + k3 = j + k * order; + + sum += basis[k2] * basis[k3]; + } + + k2 = i + j * order; + k3 = j + i * order; + alpha[k2] = sum; + alpha[k3] = sum; + + } + } + + /* Calculate beta */ + + for(i = 0; i < order; i++) { + sum = 0; + for(j = 0; j < n; j++) { + k2 = i + j * order; + sum += data[j] * basis[k2]; + } + beta[i] = sum; + } + + /* Solve for everything */ + + solve(alpha, beta, order); + + /* Now calculate and return the mean squared error */ + + sum = 0; + for(i = 0; i < n; i++) { + for(j = 0, xi = 1, val = 0; j < order; j++, xi *= x[i]) + val += beta[j] * xi; + error = val - data[i]; + sum += error * error; + } + + return sum/n; + } + + /* solve: */ + /* Uses Gaussian elimination to finish up the job for polyfit */ + + void solve(double a[], float b[], int n) { + double mag, mag2, temp; + int pivot; + int i, j, k; + + for(i = 0; i < n; i++) { + mag = 0; + pivot = -1; + + for(j = i; j < n; j++) { + mag2 = fabs(a[i + j * n]); + if(mag2 > mag) { + mag = mag2; + pivot = j; + } + } + + if(pivot == -1 || mag == 0) return; + + if(pivot != i) { + for(j = i; j < n; j++) { + temp = a[j + i * n]; + a[j + i * n] = a[j + pivot * n]; + a[j + pivot * n] = temp; + } + + temp = b[i]; + b[i] = b[pivot]; + b[pivot] = temp; + } + + mag = a[i + i * n]; + for(j = i; j < n; j++) + a[j + i * n] /= mag; + b[i] /= mag; + + for(j = 0; j < n; j++) { + if(j == i) continue; + + mag2 = a[i + j * n]; + + for(k = i; k < n; k++) + a[k + j * n] -= mag2 * a[k + i * n]; + + b[j] -= mag2 * b[i]; + } + + } + } *** src/physio/polyfit.h Wed Mar 24 17:25:21 1999 --- /SLOTH-U2/welling/fiasco/development/Fiasco4.2/src/physio//polyfit.h Tue Feb 9 12:21:09 1999 *************** *** 0 **** --- 1,10 ---- + /* Polyfit: */ + /* Arguments - */ + /* data : an array of doubles containing data points */ + /* n : # of data points */ + /* beta : array of doubles for holding the cofficients */ + /* order : order of the polynomial to fit to */ + /* Returns - */ + /* mean squared error of the data with the fitted polynomial */ + + float polyfit(float *, int *, int, float *, int); *** src/physio/resppeaks.c Wed Mar 24 17:25:21 1999 --- /SLOTH-U2/welling/fiasco/development/Fiasco4.2/src/physio//resppeaks.c Tue Feb 9 14:37:19 1999 *************** *** 0 **** --- 1,287 ---- + /************************************************************ + * * + * resppeaks.c * + * * + * Permission is hereby granted to any individual or * + * institution for use, copying, or redistribution of * + * this code and associated documentation, provided * + * that such code and documentation are not sold for * + * profit and the following copyright notice is retained * + * in the code and documentation: * + * Copyright (c) 1997 Department of Statistics, * + * Carnegie Mellon University * + * * + * Original programming by Tom Bagby, 12/98 * + ************************************************************/ + /************************************************************* + Description: + + resppeaks takes a pgh MRI dataset containing a respiratory chunk + and creates a new MRI file containing the chunk with a ramping + function cooresponding to each break in the original file. A + second chunk is also added containing magnitudes for each breath. + + Parameters: + + resppeaks [-mc mchunk] [-c chunk] [-w window_size] [-t threshold] + infile outfile + + -mc mchunk specifies the name of the new magnitude chunk + -c chunk specifies the name of the main input/output chunk + to use + -w window_size specifies the size of the window to use for + filtering the data when looking for peaks + -t threshold specifices the threshold to use for deciding + infile specifies the input dataset + output specifies the output dataset + + **********************************************************/ + + #include + #include + #include + #include "mri.h" + #include "fmri.h" + #include "misc.h" + #include "stdcrg.h" + #include "polyfit.h" + + /* Empirically derived constants */ + #define DEFAULT_WINDOW 70 + #define DEFAULT_THRESHOLD 5.0 + + /* Global parameters */ + int window; + float threshold; + + int find_next_peak(float *, int, int); + void process_data(float *, float *, float *, int); + + /* process_data: */ + /* does all the actual data processing taking the cardio data in the */ + /* input buffer and making our cooresponding ramping function in the */ + /* output buffer where n is the size of the buffers. */ + + void process_data(float *in_buffer, + float *out_buffer, + float *mag_buffer, + int n) { + + float coeffs[3], flat_buffer[4], slope, yinter, val; + int i, j, k, w2 = window/2; + int x[4]; + + float *temp_buffer = (float *)malloc(n * sizeof(float)); + + /* Set the temp and output buffer contents to 0 */ + for(i = 0; i < n; i++) { + out_buffer[i] = 0; + temp_buffer[i] = 0; + } + + /* Make a first pass over the data fitting parabolas */ + for(i = w2; i < n - w2; i++) { + polyfitlin(in_buffer + i - w2, window, coeffs, 3); + temp_buffer[i] = coeffs[2] * 200; + } + + /* Now go through and find the peaks from the fitting stage */ + /* at each peak, we seek a local minimum in the orginal */ + /* data around it. We also check to see whether we are at */ + /* a point where the data flattens out entirely. If it is */ + /* we place the peak at the middle of the flat section. */ + + for(i = w2; i < n - w2; i++) + if(temp_buffer[i-1] < temp_buffer[i] && + temp_buffer[i+1] < temp_buffer[i] && + temp_buffer[i] > threshold) { + + /* check whether we are at a flat area */ + for(j = i + 1; j < n - w2; j++) + if(fabs(in_buffer[j] - in_buffer[i]) > 15) + break; + + if(j > i + window/5) { + flat_buffer[0] = in_buffer[i - 2]; + flat_buffer[1] = in_buffer[i - 1]; + flat_buffer[2] = in_buffer[j + 1]; + flat_buffer[3] = in_buffer[j + 2]; + x[0] = i - 2; x[1] = i - 1; + x[2] = j + 1; x[3] = j + 2; + polyfit(flat_buffer, x, 4, coeffs, 3); + val = (i + j)/2.0; + + out_buffer[(i + j)/2] = 1; + mag_buffer[(i + j)/2] = fabs(coeffs[0] + + coeffs[1] * val + + coeffs[2] * val * val); + i = j; + continue; + } + + j = i; + while((in_buffer[j] > in_buffer[j+1] || + in_buffer[j] > in_buffer[j-1]) && + (j > 1 && j < n-1)) + if(in_buffer[j+1] < in_buffer[j]) j = j+1; + else j = j-1; + out_buffer[j] = 1; + mag_buffer[j] = fabs(in_buffer[j]); + } + + /* Now linearize everything */ + + i = -1; + while((i = find_next_peak(out_buffer, i, n)) != -1) { + j = find_next_peak(out_buffer, i, n); + if(j == -1) break; + slope = 1.0/(j - i - 1); + for(k = i + 1; k < j; k++) { + out_buffer[k] = (k - i - 1) * slope; + mag_buffer[k] = mag_buffer[j]; + } + + } + + /* Now we have to take care of the sections at the beginning and + end of the data that just didn't fit into our window for the + we copy the slope from the first and last sections where we were + able to identify peaks, and fill in the gaps. */ + + i = find_next_peak(out_buffer, -1, n); + j = find_next_peak(out_buffer, i, n); + for(k = 0; k < i; k++) { + out_buffer[k] = out_buffer[j - i + k]; + mag_buffer[k] = mag_buffer[i]; + } + + for(i = n - 1; i >= 0 && out_buffer[i] == 0; i--); + for(j = i - 1; j >= 0 && out_buffer[j] != 1; j--); + for(k = i + 1; k < n; k++) { + out_buffer[k] = out_buffer[j - i + k]; + mag_buffer[k] = mag_buffer[i]; + } + + } + + /* find_next_peak: */ + /* given the location of the previous peak, finds the next */ + /* peak. if none is to be found, returns -1 */ + + int find_next_peak(float *buffer, int start, int size) { + int i = start + 1; + while(i < size && buffer[i] != 1) i++; + if(i == size) return -1; + return i; + } + + int main(int argc, char *argv[]) { + MRI_Dataset *input, *output; + int chunk_size, num_elements, i, j, k; + float *input_buffer, *output_buffer, *mag_buffer; + char chunk[512]; + char mchunk[512]; + char infile[512]; + char outfile[512]; + char tag[600]; + + /* Check to see if help was requested */ + if( ( argc > 1 ) && !strcmp( argv[1], "-help" ) ) + { + if( argc == 2 ) + Help( "selecttopic" ); + else + Help( argv[2] ); + } + + /* First, parse the command line */ + + cl_scan(argc, argv); + + /* Then, get the options, input and output filenames */ + + cl_get("mc", "%option %s[%]", "magnitude", mchunk); + cl_get("c", "%option %s[%]", "physio", chunk); + cl_get("w", "%option %d[%]", DEFAULT_WINDOW, &window); + cl_get("t", "%option %f[%]", DEFAULT_THRESHOLD, &threshold); + + if(!cl_get("", "%s", infile)) { + fprintf(stderr, "%s: Input file name not given.\n", argv[0]); + exit(-1); + } + if(!cl_get("", "%s", outfile)) { + fprintf(stderr, "%s: Output file name not given.\n", argv[0]); + exit(-1); + } + + if(cl_cleanup_check()) { + int i; + fprintf(stderr, "%s: invalid argument in command line:\n", argv[0]); + for(i = 0; i < argc; i++) fprintf(stderr, "%s ", argv[i]); + fprintf(stderr, "\n"); + Help("usage"); + exit(-1); + } + + /* Now that we are done the command-line parsing, lets copy the + whole dataset and open everything up so we can start messing + around with it. If the named chunk isn't in the dataset, + then we abort. We also abort if the input and output filenames + are the same */ + + if(!strcmp(infile, outfile)) { + fprintf(stderr, "%s: input and output files must be distinct.\n", argv[0]); + exit(-1); + } + input = mri_open_dataset(infile, MRI_READ); + if(!mri_has(input, chunk)) { + fprintf(stderr, "%s: requested chunk %s not present.\n", argv[0], chunk); + exit(-1); + } + strcpy(tag, chunk); + strcat(tag, ".dimensions"); + if (strcmp(mri_get_string(input,tag), "t")) { + fprintf(stderr,"%s: input chunk does not have dimensions 't'\n", argv[0]); + exit(-1); + } + output = mri_copy_dataset(outfile, input); + + /* Get the length of the chunk */ + strcpy(tag, chunk); + strcat(tag, ".extent.t"); + chunk_size = mri_get_int(input, tag); + + /* Now we read everything in and get the output buffer. Since the physio + data sets should not be ridiculously huge, everything gets read in + at the same time. */ + + mri_create_chunk(output, "magnitude"); + + strcpy(tag, mchunk); + strcat(tag, ".datatype"); + mri_set_string(output, tag, "float32"); + + strcpy(tag, mchunk); + strcat(tag, ".dimensions"); + mri_set_string(output, tag, "t"); + + strcpy(tag, mchunk); + strcat(tag, ".extent.t"); + mri_set_int(output, tag, chunk_size); + + input_buffer = mri_get_chunk(input, chunk, chunk_size, 0, MRI_FLOAT); + output_buffer = mri_get_chunk(output, chunk, chunk_size, 0, MRI_FLOAT); + mag_buffer = mri_get_chunk(output, "magnitude", chunk_size, 0, MRI_FLOAT); + + /* Then we process the data, write out the results and close our + files. Then there is much happiness and joy. */ + + process_data(input_buffer, output_buffer, mag_buffer, chunk_size); + mri_set_chunk(output, chunk, chunk_size, 0, MRI_FLOAT, output_buffer); + mri_set_chunk(output, "magnitude", chunk_size, 0, MRI_FLOAT, mag_buffer); + + mri_close_dataset(input); + mri_close_dataset(output); + + exit(0); + } *** src/physio/resppeaks_help.help Wed Mar 24 17:25:21 1999 --- /SLOTH-U2/welling/fiasco/development/Fiasco4.2/src/physio//resppeaks_help.help Tue Feb 9 14:44:26 1999 *************** *** 0 **** --- 1,31 ---- + *Usage + + Description: + + resppeaks takes a pgh MRI dataset containing a respiratory chunk + and creates a new MRI file containing the chunk with a ramping + function cooresponding to each break in the original file. + Specifically, the output chunk values run linearly from 0.0 + immediately following one inhalation maximum to 1.0 the moment + before the next. Estimates are made for the values before and + after the first and last peaks. A second chunk (named "magnitude") + is also added containing magnitudes for each breath. + + The input chunk must be a pure time series, that is, it's + dimensions must be "t". + + Parameters: + + resppeaks [-mc mchunk] [-c chunk] [-w window_size] [-t threshold] + infile outfile + + -mc mchunk specifies the name of the new magnitude chunk + -c chunk specifies the name of the main input/output chunk + to use + -w window_size specifies the size of the window to use for + filtering the data when looking for peaks + -t threshold specifices the threshold to use for deciding + infile specifies the input dataset + output specifies the output dataset + + Index:src/brain/generate_vmpfx_in.c 35,36c35 < -fixed names of factor levels to hold fixed, as a comma-delimed list < (default "NA") --- > -fixed names of factor levels to hold fixed; this is required 38a38 > -null flag indicating that only the null model is to be run 79c79 < static char rcsid[] = "$Id: generate_vmpfx_in.c,v 1.8 1998/09/30 22:28:45 welling Exp $"; --- > static char rcsid[] = "$Id: generate_vmpfx_in.c,v 1.9 1999/03/18 23:43:41 welling Exp $"; 89a90 > static char condition_map_string[STRING_LENGTH]; 100a102,103 > static int some_images_NA= 0; > static int null_model_only= 0; 198,199c201,202 < char* this_cond= strtok(fixed_cond_str," ,"); < char* runner= fixed_cond_out_string; --- > char* this_cond; > 201c204,208 < *runner= 0; --- > /* The missing condition is always fixed by assertion */ > cond_def_table[0]->fixed= 1; > n_fixed_conditions++; > > this_cond= strtok(fixed_cond_str," ,"); 213,217d219 < if ((int)(runner-fixed_cond_out_string) > STRING_LENGTH-64) < Abort("%s: fixed cond string overflow!\n",progname); < if (runner != fixed_cond_out_string) *runner++= ','; < sprintf(runner,"%d",i); < runner += strlen(runner); 321a324,430 > static int test_and_filter_conditions() > { > int* cond_guesses; > int i; > int icond; > char* runner= fixed_cond_out_string; > > /* Generate consistent guesses about the conditions in effect > * for each image. > */ > if (!(cond_guesses= (int*)malloc(nimages*sizeof(int)))) > Abort("%s: unable to allocate %d bytes!\n",progname,nimages*sizeof(int)); > for (i=0; i for (icond=0; icond if (split_table[icond].cond != 0) { > int image= split_table[icond].image; > > if (cond_guesses[image]==0) cond_guesses[image]= split_table[icond].cond; > else { > /* The following test determines if we have already seen > * a slice in this image with a non-NA cond. This would > * violate the requirements for vmpfx which we choose to > * enforce. > */ > if (cond_guesses[image] != split_table[icond].cond) return 0; > } > } > } > > /* Go through and refill any slices marked NA (condition 0) > * with the value appropriate for the rest of that image. > * This allows some NAs to persist, but for the whole image > * only. > */ > for (icond=0; icond if (split_table[icond].cond==0) > split_table[icond].cond= cond_guesses[ split_table[icond].image ]; > } > > /* Make a note of whether or not any trace of the NA > * condition remains. > */ > some_images_NA= 0; > for (i=0; i if (cond_guesses[i]==0) some_images_NA= 1; > > /* Make sure we have enough fixed conditions to get by */ > if ((some_images_NA && n_fixed_conditions<1) > || (!some_images_NA && n_fixed_conditions<2)) > Abort("%s: at least one fixed condition must exist in experiment!\n", > progname); > > /* Build the output list of fixed conditions, including > * or excluding the NA condition as appropriate. > */ > runner= fixed_cond_out_string; > *runner= 0; /* null terminate */ > if (cond_def_table[0]->fixed && some_images_NA) { > sprintf(runner,"0 "); > runner += 2; > } > for (i=1; i if (cond_def_table[i]->fixed) { > if ((int)(runner-fixed_cond_out_string) > STRING_LENGTH-64) > Abort("%s: fixed cond string overflow!\n",progname); > if (runner != fixed_cond_out_string) *runner++= ' '; > sprintf(runner,"%d", some_images_NA ? i : i-1); > runner += strlen(runner); > } > } > > /* Build a list of Fiasco condition numbers corresponding to > * the valid range of Brain condition numbers. Format is > * a space-separated list of f/(brainCondNumber).(fiascoCondNumber) > * triples, where f is replaced by v for non-fixed conditions. > * It needn't include all fiascoCondNumbers, since the NA > * condition may be dropped. The first brainCondNumber > * will always be zero. > */ > runner= condition_map_string; > *runner= 0; /* null terminate */ > if (some_images_NA) { > for (i=0; i if ((int)(runner-condition_map_string) > STRING_LENGTH-64) > Abort("%s: map cond string overflow!\n",progname); > sprintf(runner,"%c/%d.%d ", > (cond_def_table[i]->fixed ? 'f' : 'v'), > i,i); > runner += strlen(runner); > } > } > else { > for (i=1; i if ((int)(runner-condition_map_string) > STRING_LENGTH-64) > Abort("%s: map cond string overflow!\n",progname); > sprintf(runner,"%c/%d.%d ", > (cond_def_table[i]->fixed ? 'f' : 'v'), > i-1,i); > runner += strlen(runner); > } > } > > free(cond_guesses); > > return 1; > } > 325c434,435 < float current_start= split_table[0].start_time; --- > int cond_out; > int current_start; 327a438 > current_start= 0; 330,331c441,445 < printf("%d %f %f\n", current_condition, current_start, < split_table[i].start_time - current_start); --- > cond_out= (some_images_NA) ? > current_condition : current_condition-1; > printf("%d %d %d\n", cond_out, > current_start, split_table[i].image-current_start); > current_start= split_table[i].image; 333d446 < current_start= split_table[i].start_time; 337,338c450,451 < printf("%d %f %f\n", current_condition, current_start, < split_table[(nimages*nslices)-1].start_time - current_start); --- > printf("%d %d %d\n", cond_out, current_start, > split_table[(nimages*nslices)-1].image + 1 - current_start); 344a458,465 > > /* Check for conditions forbidden with the -null flag */ > if (null_model_only && > (strstr(instring,"****CONDITION_TABLE****") > || strstr(instring,"****NUM_CONDITIONS****") > || strstr(instring,"****FIXED_CONDITIONS****"))) > Abort("%s: prototype file is inappropriate for -null flag\n",progname); > 354c475,483 < sprintf(tstring,"%d",n_conditions); --- > > sprintf(tstring,"%d",nimages); > (void)check_and_substitute(instring,"****NUM_IMAGES****",tstring); > > sprintf(tstring,"%d",nslices); > (void)check_and_substitute(instring,"****NUM_SLICES****",tstring); > > if (some_images_NA) sprintf(tstring,"%d",n_conditions); > else sprintf(tstring,"%d",n_conditions-1); 355a485 > 357a488,489 > (void)check_and_substitute(instring,"****CONDITION_MAP****", > condition_map_string); 412,414c544,552 < cl_get( "split", "%option %s[%]", "split", split_fname); < cl_get( "cond", "%option %s[%s]", "conditions", cond_fname); < cl_get( "fixed", "%option %s[%s]", "NA", fixed_conditions); --- > null_model_only= cl_present("null"); > if (!null_model_only) { > cl_get( "split", "%option %s[%]", "split", split_fname); > cl_get( "cond", "%option %s[%s]", "conditions", cond_fname); > if (!cl_get( "fixed", "%option %s", fixed_conditions)) { > fprintf(stderr,"%s: -fixed option required.\n",argv[0]); > Help( "usage" ); > } > } 429,432c567 < /* Parse condition file */ < if (!(cond_file= fopen(cond_fname,"r"))) < Abort("%s: unable to open condition file <%s>!\n", < argv[0], cond_fname); --- > if (!null_model_only) { 434,460c569,604 < if (!process_cond_file(cond_file)) < Abort("%s: error processing conditions file <%s>!\n", < argv[0],cond_fname); < < if (n_conditions>MAX_SAFE_CONDITIONS) < fprintf(stderr, < "%s: *****WARNING***** vmpfx sometimes fails to converge above %d conditions; you have %d!\n", < argv[0], MAX_SAFE_CONDITIONS, n_conditions); < < if (fclose(cond_file)) < Error("%s: failed to close conditions file!\n",argv[0]); < < /* Note fixed conditions ***/ < if (!parse_fixed_cond(fixed_conditions)) < Abort("%s: failed to parse fixed condition list\n", < progname); < < /* Parse split file */ < if (!(split_file= fopen(split_fname,"r"))) < Abort("%s: unable to open split file <%s>!\n", < argv[0],split_fname); < < if (!process_split_file(split_file)) < Abort("%s: error processing split file <%s>!\n",argv[0],split_fname); < < if (fclose(split_file)) < Error("%s: failed to close split file!\n",argv[0]); --- > /* Parse condition file */ > if (!(cond_file= fopen(cond_fname,"r"))) > Abort("%s: unable to open condition file <%s>!\n", > argv[0], cond_fname); > > if (!process_cond_file(cond_file)) > Abort("%s: error processing conditions file <%s>!\n", > argv[0],cond_fname); > > if (n_conditions>MAX_SAFE_CONDITIONS) > fprintf(stderr, > "%s: *****WARNING***** vmpfx sometimes fails to converge above %d conditions; you have %d!\n", > argv[0], MAX_SAFE_CONDITIONS, n_conditions); > > if (fclose(cond_file)) > Error("%s: failed to close conditions file!\n",argv[0]); > > /* Note fixed conditions ***/ > if (!parse_fixed_cond(fixed_conditions)) > Abort("%s: failed to parse fixed condition list\n", > progname); > > /* Parse split file */ > if (!(split_file= fopen(split_fname,"r"))) > Abort("%s: unable to open split file <%s>!\n",argv[0],split_fname); > > if (!process_split_file(split_file)) > Abort("%s: error processing split file <%s>!\n",argv[0],split_fname); > > if (fclose(split_file)) > Error("%s: failed to close split file!\n",argv[0]); > > if (!test_and_filter_conditions()) > Abort("%s: this experimental design is not compatible with vmpfx!\n", > argv[0]); > } Index:src/brain/generate_vmpfx_in_help.help 20a21,23 > -iai inter-acquisition interval value (default 3.0) > -null specifies that the given prototype file requests only the > null model, and thus does not need split, cond, or fixed info 23d25 < -iai inter-acquisition interval value (default 3.0) 25c27,28 < (default "NA") --- > (required unless -null is specified) > Index:src/brain/vmpfx.c 987c987 < # if defined(NON_ADDITIVE_PASTE) --- > #if defined(NON_ADDITIVE_PASTE) 995c995 < # endif --- > #endif 1776c1776 < # endif --- > #endif 1804c1804 < #if defined( USE_SIGNALS ) --- > #if defined( USE_SIGNALS ) 1806c1806 < #endif --- > #endif 1923c1923 < #if defined(MONITOR) && (MONITOR > 0) --- > #if defined(MONITOR) && (MONITOR > 0) 1931c1931 < #endif --- > #endif 1933c1933 < #if DIAGNOSTIC >= 5 --- > #if DIAGNOSTIC >= 5 1935c1935 < #endif --- > #endif 2016c2016 < #if DIAGNOSTIC >= 7 --- > #if DIAGNOSTIC >= 7 2023c2023 < #endif --- > #endif 2074c2074 < #if DIAGNOSTIC >= 6 --- > #if DIAGNOSTIC >= 6 2076c2076 < #endif --- > #endif 2080c2080 < #if DIAGNOSTIC >= 8 --- > #if DIAGNOSTIC >= 8 2087c2087 < #endif --- > #endif 2162c2162 < #if DIAGNOSTIC >= 9 --- > #if DIAGNOSTIC >= 9 2164c2164 < #endif --- > #endif 2173c2173 < #if DIAGNOSTIC >= 1 --- > #if DIAGNOSTIC >= 7 2175c2175 < #endif --- > #endif 2542c2542 < #if DIAGNOSTIC >= 20 --- > #if DIAGNOSTIC >= 20 2567c2567 < #endif --- > #endif 2759c2759 < #if DIAGNOSTIC >= 20 --- > #if DIAGNOSTIC >= 20 2777c2777 < #endif --- > #endif 2846c2846 < #if DIAGNOSTIC >= 20 --- > #if DIAGNOSTIC >= 20 2855c2855 < #if DIAGNOSTIC >= 20 --- > #if DIAGNOSTIC >= 20 2885c2885 < #endif --- > #endif 2897c2897 < #endif --- > #endif 3147c3147 < #if !defined( FIX_SIGMA ) --- > #if !defined( FIX_SIGMA ) 3156c3156 < #else --- > #else 3158c3158 < #define SIGMA_REFINE 3 /* ATTN: SIGMA_REFINE is TEMPORARY and PROVISIONAL */ --- > #define SIGMA_REFINE 3 /* ATTN: SIGMA_REFINE is TEMPORARY and PROVISIONAL */ 3180c3180 < #endif --- > #endif 3201c3201 < #if DIAGNOSTIC >= 20 --- > #if DIAGNOSTIC >= 20 3217c3217 < #if DIAGNOSTIC >= 20 --- > #if DIAGNOSTIC >= 20 3253c3253 < #endif --- > #endif 3257c3257 < #endif --- > #endif 3371c3371 < #if DIAGNOSTIC>=2 --- > #if DIAGNOSTIC>=2 3374c3374 < #endif --- > #endif 3380c3380 < #if DIAGNOSTIC>=2 --- > #if DIAGNOSTIC>=2 3382c3382 < #endif --- > #endif 3386c3386 < #if DIAGNOSTIC >= 1 --- > #if DIAGNOSTIC >= 1 3389c3389 < #endif --- > #endif 3693c3693 < #if DIAGNOSTIC >= 10 --- > #if DIAGNOSTIC >= 10 3701c3701 < #endif --- > #endif 3827c3827 < #if DIAGNOSTIC >= 10 --- > #if DIAGNOSTIC >= 10 3851c3851 < #endif --- > #endif 6226c6226 < # if defined(USE_MRIFUNC) --- > #if defined(USE_MRIFUNC) 6250c6250 < # else --- > #else 6254c6254 < # endif --- > #endif Index:src/brain/vmpfx_proto.t 24c24 < design.images 0 --- > design.images 1 28a29,30 > #FiascoCondMap ****CONDITION_MAP**** > 53,54c55 < #prior.sigmasq 1.6 100.0 < prior.sigmasq 5.0 1.5 --- > prior.sigmasq 1.8 200.0 Index:src/csh/FIASCO 22c22 < echo '$Id: FIASCO,v 1.12 1998/12/03 00:30:19 welling Exp $'\ --- > echo '$Id: FIASCO,v 1.13 1999/03/24 22:16:20 welling Exp $'\ 29c29 < setenv FIASCO_PATCHLVL 1 --- > setenv FIASCO_PATCHLVL 2 Index:src/csh/cg_brain.csh 19c19 < echo '#$Id: cg_brain.csh,v 1.7 1998/09/24 22:23:25 welling Exp $' --- > echo '#$Id: cg_brain.csh,v 1.8 1999/03/24 21:37:37 welling Exp $' 28a29,30 > smooth_missing -i $3.mri -h ${3}_s.mri > rm $3.mri $3.dat 32c34 < $F_NSLICE -droot ${3:t} -dpath ${3:h}/ -split in/$F_SPLIT_NEW \ --- > $F_NSLICE -droot ${3:t}_s -dpath ${3:h}/ -split in/$F_SPLIT_NEW \ 36a39,45 > > # Clean up some storage, and recreate the original (unsmoothed) > # intermediate file. We do this for historical compatibility. > rm ${3}_s.mri ${3}_s.dat > echo '#'`date` > permute -memlimit 32000000 -input $1.mri -headerout $3.mri \ > -dataout .dat -chunk images -order vtxyz Index:src/csh/cg_brain_stats.csh 22c22 < echo '#$Id: cg_brain_stats.csh,v 1.8 1998/09/29 16:57:45 welling Exp $' --- > echo '#$Id: cg_brain_stats.csh,v 1.10 1999/03/24 21:39:13 welling Exp $' 30c30 < echo 'Expected vmpfx output does not exist!' --- > echo 'Expected vmpfx output does not exist.' 80,81c80,81 < # conditions that are not fixed. This requires manipulating the < # fixed condition and condition table lists. --- > # conditions that are not fixed. This is done using special > # translation info saved earlier by generate_vmpfx_in. 83,96c83 < # need a local variable to hold fixed cond info < set vmpfx_fixed = ( $F_VMPFX_FIXED ) < < # make a list of ints 0 <= n < n_resp < set ncond_list = 0 < set count = 1 < while ($count < $n_resp) < set ncond_list = ( $ncond_list $count ) < @ count = $count + 1 < end < < # extract stats for experimental conditions that do not appear in < # the fixed condition list. We build up a list of non-fixed conditions < # at the same time; the first bit makes unfixed_list an empty variable. --- > # Set up a list of unfixed Fiasco indices 99d85 < foreach cond ($ncond_list) 101,118c87,100 < set cond_string = `cat $F_SPLIT_COND | grep -E '^[ ]*'$cond` < set cond_fixed = 0 < set cond_index = 1 < while ($cond_index <= $#cond_string) < set fixed_index = 1 < while ($fixed_index <= $#vmpfx_fixed) < if ($vmpfx_fixed[$fixed_index] == $cond_string[$cond_index]) then < set cond_fixed = 1 < break < endif < @ fixed_index = $fixed_index + 1 < end < if ( $cond_fixed ) break < @ cond_index = $cond_index + 1 < end < if ( ! $cond_fixed ) then < set unfixed_list = ( $unfixed_list $cond ) < @ cond_offset = $n_baselevel + $n_drift + $n_shape + $n_noise + $cond + 1 --- > # need the map between brain and Fiasco conds, generated > # by generate_vmpfx_in. > set cond_map = `grep FiascoCondMap $F_VMPFX_INPUT` > shift cond_map > > # Walk the condition map, generating stats for non-fixed > # conditions. > foreach mapentry ($cond_map) > if (${mapentry:h} == v) then > set pair = ${mapentry:t} > set bindex = ${pair:r} > set findex = ${pair:e} > set unfixed_list = ( $unfixed_list $findex ) > @ cond_offset = $n_baselevel + $n_drift + $n_shape + $n_noise + $bindex + 1 121c103 < $F_STAT_MPATH/Resp.$cond.Mean --- > $F_STAT_MPATH/Resp.$findex.Mean 123,127c105,107 < $F_STAT_SPATH/Resp.$cond.Stdv < # mri_rpn_math '$1,$2,/' -o $F_STAT_TPATH/Resp.$cond.Tmap \ < # $F_STAT_MPATH/Resp.$cond.Mean $F_STAT_SPATH/Resp.$cond.Stdv < mri_rpn_math '$1,$2,*' -o $F_STAT_TPATH/Bayes.PResp.$cond \ < $F_STAT_MPATH/Resp.$cond.Mean $F_STAT_SPATH/Bayes.P --- > $F_STAT_SPATH/Resp.$findex.Stdv > mri_rpn_math '$1,$2,*' -o $F_STAT_TPATH/Bayes.PResp.$findex \ > $F_STAT_MPATH/Resp.$findex.Mean $F_STAT_SPATH/Bayes.P Index:src/csh/epi.defaults.csh 17c17 < echo '#$Id: epi.defaults.csh,v 1.19 1998/09/23 22:13:08 welling Exp $' --- > echo '#$Id: epi.defaults.csh,v 1.20 1999/03/24 21:34:19 welling Exp $' 274a275,283 > # > # > #spline_detrend.csh > # > setenv F_VMPFX_NPROTO ${FIASCO}/../../src/brain/vmpfx_null_proto.t > # > setenv F_VMPFX_NINPUT in/vmpfx_null_in.t > # image acquisition interval in seconds > setenv F_IAI 3.0 Index:src/csh/epi.local.csh 17c17 < # $Id: epi.local.csh,v 1.8 1998/09/23 21:47:10 welling Exp $ --- > # $Id: epi.local.csh,v 1.9 1999/03/24 21:35:40 welling Exp $ 49a50,53 > # > # image acquisition interval in seconds > setenv F_IAI 3.0 > # Index:src/csh/estimate_z_motion.csh 24c24 < echo '# $Id: estimate_z_motion.csh,v 1.2 1998/08/28 20:18:07 welling Exp $' --- > echo '# $Id: estimate_z_motion.csh,v 1.3 1999/02/11 21:52:38 welling Exp $' 26c26,28 < # Pull in some parameters for smregpar --- > # Pull in some parameters for smregpar. This will smash in/split, so > # save and restore it. > cp in/split in/e_z_motion_split_sav 27a30,31 > rm in/split > mv in/e_z_motion_split_sav in/split Index:src/csh/spiral.defaults.csh 17c17 < echo '$Id: spiral.defaults.csh,v 1.11 1998/09/23 22:23:09 welling Exp $' --- > echo '$Id: spiral.defaults.csh,v 1.12 1999/03/24 21:34:32 welling Exp $' 145a146,154 > # > # > #spline_detrend.csh > # > setenv F_VMPFX_NPROTO ${FIASCO}/../../src/brain/vmpfx_null_proto.t > # > setenv F_VMPFX_NINPUT in/vmpfx_null_in.t > # image acquisition interval in seconds > setenv F_IAI 3.0 Index:src/csh/spiral.local.csh 17c17 < #$Id: spiral.local.csh,v 1.7 1998/09/23 22:28:12 welling Exp $ --- > #$Id: spiral.local.csh,v 1.8 1999/03/24 21:36:10 welling Exp $ 41a42,45 > # > # image acquisition interval in seconds > setenv F_IAI 3.0 > # Index:src/csh/ts.defaults.csh 17c17 < echo '#$Id: ts.defaults.csh,v 1.19 1998/09/23 22:12:51 welling Exp $' --- > echo '#$Id: ts.defaults.csh,v 1.20 1999/03/24 21:34:48 welling Exp $' 308a309,317 > # > # > #spline_detrend.csh > # > setenv F_VMPFX_NPROTO ${FIASCO}/../../src/brain/vmpfx_null_proto.t > # > setenv F_VMPFX_NINPUT in/vmpfx_null_in.t > # image acquisition interval in seconds > setenv F_IAI 3.0 Index:src/csh/ts.local.csh 17c17 < # $Id: ts.local.csh,v 1.9 1998/09/23 21:48:34 welling Exp $ --- > # $Id: ts.local.csh,v 1.10 1999/03/24 21:35:55 welling Exp $ 52a53,56 > # > # image acquisition interval in seconds > setenv F_IAI 3.0 > # Index:src/estireg/estireg_help.help 86,87c86,87 < "Improved Image Registration Using Fourier Interpolation", to < appear: Magnetic Resonance in Medicine. --- > "Improved Image Registration Using Fourier Interpolation", > Magnetic Resonance in Medicine 36-6 (Dec 1996) Index:src/fmri/fft2d.c 50c50 < static char rcsid[] = "$Id: fft2d.c,v 1.7 1998/09/25 16:03:31 welling Exp $"; --- > static char rcsid[] = "$Id: fft2d.c,v 1.9 1998/12/01 19:40:41 welling Exp $"; 131c131 < static int check_plan_x( int xdim, int ydim, fftw_direction xdir ) --- > static void check_plan_x( int xdim, int ydim, fftw_direction xdir ) 161c161 < static int check_plan_y( int xdim, int ydim, fftw_direction ydir ) --- > static void check_plan_y( int xdim, int ydim, fftw_direction ydir ) 191c191 < static int check_plan_2d( int xdim, int ydim, fftw_direction dir ) --- > static void check_plan_2d( int xdim, int ydim, fftw_direction dir ) 221c221 < void data_to_temp_2d( FComplex** data, long nx, long ny ) --- > static void data_to_temp_2d( FComplex** data, long nx, long ny ) 264c264 < void data_from_temp_2d( FComplex** data, long nx, long ny, double scale ) --- > static void data_from_temp_2d( FComplex** data, long nx, long ny, double scale ) 307c307 < void data_to_temp_x( FComplex** data, long nx, long ny ) --- > static void data_to_temp_x( FComplex** data, long nx, long ny ) 330c330 < void data_from_temp_x( FComplex** data, long nx, long ny, double scale ) --- > static void data_from_temp_x( FComplex** data, long nx, long ny, double scale ) 353c353 < void data_to_temp_y( FComplex** data, long nx, long ny ) --- > static void data_to_temp_y( FComplex** data, long nx, long ny ) 376c376 < void data_from_temp_y( FComplex** data, long nx, long ny, double scale ) --- > static void data_from_temp_y( FComplex** data, long nx, long ny, double scale ) Index:src/fmri/fft3d.c 55c55 < static char rcsid[] = "$Id: fft3d.c,v 1.3 1998/10/01 22:08:33 welling Exp $"; --- > static char rcsid[] = "$Id: fft3d.c,v 1.5 1998/12/01 19:41:28 welling Exp $"; 147c147 < static int check_plan_x( int xdim, int ydim, int zdim, fftw_direction xdir ) --- > static void check_plan_x( int xdim, int ydim, int zdim, fftw_direction xdir ) 178c178 < static int check_plan_y( int xdim, int ydim, int zdim, fftw_direction ydir ) --- > static void check_plan_y( int xdim, int ydim, int zdim, fftw_direction ydir ) 211c211 < static int check_plan_z( int xdim, int ydim, int zdim, fftw_direction zdir ) --- > static void check_plan_z( int xdim, int ydim, int zdim, fftw_direction zdir ) 244c244 < static int check_plan_xy( int xdim, int ydim, int zdim, fftw_direction dir ) --- > static void check_plan_xy( int xdim, int ydim, int zdim, fftw_direction dir ) 277c277 < static int check_plan_yz( int xdim, int ydim, int zdim, fftw_direction dir ) --- > static void check_plan_yz( int xdim, int ydim, int zdim, fftw_direction dir ) 310c310 < static int check_plan_xyz( int xdim, int ydim, int zdim, fftw_direction dir ) --- > static void check_plan_xyz( int xdim, int ydim, int zdim, fftw_direction dir ) 343c343 < void data_to_temp_x( FComplex* data, long nx, long ny, long nz ) --- > static void data_to_temp_x( FComplex* data, long nx, long ny, long nz ) 377,378c377,378 < void data_from_temp_x( FComplex* data, long nx, long ny, long nz, < double scale ) --- > static void data_from_temp_x( FComplex* data, long nx, long ny, long nz, > double scale ) 412c412 < void data_to_temp_y( FComplex* data, long nx, long ny, long nz ) --- > static void data_to_temp_y( FComplex* data, long nx, long ny, long nz ) 446,447c446,447 < void data_from_temp_y( FComplex* data, long nx, long ny, long nz, < double scale ) --- > static void data_from_temp_y( FComplex* data, long nx, long ny, long nz, > double scale ) 481c481 < void data_to_temp_z( FComplex* data, long nx, long ny, long nz ) --- > static void data_to_temp_z( FComplex* data, long nx, long ny, long nz ) 515,516c515,516 < void data_from_temp_z( FComplex* data, long nx, long ny, long nz, < double scale ) --- > static void data_from_temp_z( FComplex* data, long nx, long ny, long nz, > double scale ) 550c550 < void data_to_temp_xy( FComplex* data, long nx, long ny, long nz ) --- > static void data_to_temp_xy( FComplex* data, long nx, long ny, long nz ) 618,619c618,619 < void data_from_temp_xy( FComplex* data, long nx, long ny, long nz, < double scale ) --- > static void data_from_temp_xy( FComplex* data, long nx, long ny, long nz, > double scale ) 687c687 < void data_to_temp_yz( FComplex* data, long nx, long ny, long nz ) --- > static void data_to_temp_yz( FComplex* data, long nx, long ny, long nz ) 755,756c755,756 < void data_from_temp_yz( FComplex* data, long nx, long ny, long nz, < double scale ) --- > static void data_from_temp_yz( FComplex* data, long nx, long ny, long nz, > double scale ) 824c824 < void data_to_temp_xz( FComplex* data, long nx, long ny, long nz ) --- > static void data_to_temp_xz( FComplex* data, long nx, long ny, long nz ) 892,893c892,893 < void data_from_temp_xz( FComplex* data, long nx, long ny, long nz, < double scale ) --- > static void data_from_temp_xz( FComplex* data, long nx, long ny, long nz, > double scale ) 961c961 < void data_to_temp_xyz( FComplex* data, long nx, long ny, long nz ) --- > static void data_to_temp_xyz( FComplex* data, long nx, long ny, long nz ) 1106,1107c1106,1107 < void data_from_temp_xyz( FComplex* data, long nx, long ny, long nz, < double scale ) --- > static void data_from_temp_xyz( FComplex* data, long nx, long ny, long nz, > double scale ) Index:src/fmri/fmri.h 54a55,60 > > /* Header for quaternion.c */ > #include "quaternion.h" > > /* Header for fshrot3d.c */ > #include "fshrot3d.h" Index:src/imtops/imtops.c 88c88 < static char rcsid[] = "$Id: imtops.c,v 1.13 1998/09/11 18:51:37 welling Exp $"; --- > static char rcsid[] = "$Id: imtops.c,v 1.14 1998/12/01 19:43:48 welling Exp $"; 404c404 < for (i=0; i < imht; i++)fprintf(out,"%*0x",2,(unsigned int)((256*i)/imht)); --- > for (i=0; i < imht; i++)fprintf(out,"%02x",(unsigned int)((256*i)/imht)); 418c418 < --- > Index:src/ireg/ireg_help.help 45c45 < Default value is "ireg.dat". --- > Default value is ".dat". 75,76c75,76 < "Improved Image Registration Using Fourier Interpolation", to < appear: Magnetic Resonance in Medicine. --- > "Improved Image Registration Using Fourier Interpolation", > Magnetic Resonance in Medicine 36-6 (Dec 1996) Index:src/libcrg/buildhelp.c 172,173c172,175 < line = (char *) ecalloc( MAX_LINE + sbuf.st_size, sizeof(char) ); < --- > if (stdin_flg) > line = (char *) ecalloc( MAX_LINE, sizeof(char) ); > else > line = (char *) ecalloc( MAX_LINE + sbuf.st_size, sizeof(char) ); Index:src/libmri/libmri.c 40c40 < static char rcsid[] = "$Id: libmri.c,v 1.9 1998/09/03 20:42:32 welling Exp $"; --- > static char rcsid[] = "$Id: libmri.c,v 1.12 1999/03/04 19:33:34 welling Exp $"; 49c49 < static FILE *log = NULL; /* for debugging */ --- > static FILE *logfile = NULL; /* for debugging */ 117,118c117,118 < if (log == NULL) < log = fopen("libmri.log", "a"); --- > if (logfile == NULL) > logfile = fopen("libmri.log", "a"); 212c212 < MRI_Chunk *ch; --- > MRI_Chunk *ch, *nch; 216a217 > char chunk_filename[MRI_MAX_FILENAME_LENGTH+1]; 273c274 < for (total = 0; total < ch->size; total += BUFFER_SIZE) --- > if (ds->mode == MRI_READ) 275,279c276,307 < n = ch->size - total; < if (n > BUFFER_SIZE) < n = BUFFER_SIZE; < ptr = mri_get_chunk(ds, ch->name, n, total, MRI_RAW); < mri_set_chunk(nds, ch->name, n, total, MRI_RAW, ptr); --- > /* we make a lazy copy of the original chunk */ > for (nch = nds->chunks; nch != NULL; nch = nch->next) > if (strcmp(nch->name, ch->name) == 0) > break; > if (nch == NULL) > { > mri_report_error(ds, "libmri: internal error -- chunk name not copied correctly\n"); > abort(); > } > /* make sure we have an absolute filename (one with a '/' in it) > so that the call to GetFile does not try to relativize it to > the new dataset name */ > if (ch->file->name[0] == '/') > strcpy(chunk_filename, ch->file->name); > else > { > if (getcwd(chunk_filename, MRI_MAX_FILENAME_LENGTH) == NULL) > { > mri_report_error(ds, "libmri: cannot get current directory name\n"); > abort(); > } > strcat(chunk_filename, "/"); > strcat(chunk_filename, ch->file->name); > } > nch->actual_file = GetFile(nds, chunk_filename); > nch->actual_file->external = TRUE; > nch->actual_datatype = ch->datatype; > nch->actual_little_endian = ch->little_endian; > nch->actual_offset = ch->offset; > nch->actual_size = ch->size; > nch->modified = TRUE; > nds->recompute_positions = TRUE; 280a309,318 > else > /* we copy the chunks immediately */ > for (total = 0; total < ch->size; total += BUFFER_SIZE) > { > n = ch->size - total; > if (n > BUFFER_SIZE) > n = BUFFER_SIZE; > ptr = mri_get_chunk(ds, ch->name, n, total, MRI_RAW); > mri_set_chunk(nds, ch->name, n, total, MRI_RAW, ptr); > } 3491c3529 < fprintf(log, "%s: ", ts); --- > fprintf(logfile, "%s: ", ts); 3494c3532 < vfprintf(log, fmt, args); --- > vfprintf(logfile, fmt, args); 3496c3534 < fflush(log); --- > fflush(logfile); Index:src/mri_util/mri_subsample.c 93c93 < static char rcsid[] = "$Id: mri_subsample.c,v 1.4 1998/07/31 21:55:15 welling Exp $"; --- > static char rcsid[] = "$Id: mri_subsample.c,v 1.5 1998/10/14 18:28:45 welling Exp $"; 426c426 < ratio= ((float)selected_extent)/((float)new_extent); --- > ratio= ((double)selected_extent)/((double)new_extent); Index:src/mri_util/slave_splus.h 18a19,22 > #ifdef LINUX > #include > #endif > Index:src/smregpar/smregpar.c 43c43 < static char rcsid[] = "$Id: smregpar.c,v 1.7 1998/08/28 20:13:32 welling Exp $"; --- > static char rcsid[] = "$Id: smregpar.c,v 1.8 1999/01/05 23:34:53 welling Exp $"; 259a260,262 > > /* A needed heading */ > fprintf(stdout,"Slice Median MSE IQR\n"); 296a300 > fprintf(stdout," %d %g %g\n",z,medianmse,iqrmse); Index:src/smregpar/smregpar_help.help 5c5,6 < to register image data. --- > to register image data. In addition, a slicewise table of median > MSE and inter-quartile range values is written to standard output. Index:src/util/libmdbg.c 1a2 > #include 15c16 < static char rcsid[] = "$Id: libmdbg.c,v 1.3 1998/06/04 15:35:52 welling Exp $"; --- > static char rcsid[] = "$Id: libmdbg.c,v 1.4 1999/03/03 00:54:58 welling Exp $"; 51a53,60 > void *my_strdup( const char* s ) > { > char* sdup; > sdup= (char*)my_malloc( strlen(s)+1 ); > (void)strcpy( sdup, s ); > return sdup; > } > 55a65 > int i; 71a82,85 > /* Zero out the freed memory, in the hopes of tripping up anything > * that tries to access it later. > */ > for (i=0; isize; i++) *(p + i)= 0; Index:src/util/mdbg.h 1a2,3 > void* my_malloc(size_t size); > 2a5,6 > void my_free(void* ptr); > 3a8,9 > void* my_realloc(void* ptr, size_t size); > 4a11,14 > void* my_calloc(size_t nelem, size_t elsize); > > #define strdup my_strdup > char *my_strdup(const char *s); Index:doc/fiasco_executables.html 18,19c18,20 < a general linear model package (in glm.c) and a smoothing < package (in smoother.c). --- > a general linear model package (in glm.c), a smoothing > package (in smoother.c), and a quaternions package > (in quaternion.c). 44a46,47 >
cardiopeaks >
finds peaks in cardiac physiological data. 80a84,85 >
ireg3d >
motion correction in 3D 98a104,105 >
mri_pad >
Enlarge one dimension of a Pgh MRI dataset by padding. 145a153,154 >
resppeaks >
finds peaks and cycle depths in respiratory physiological data. 147a157,158 >
smooth_missing >
linearly interpolates over missing data in time series. 173c184 < Last modified: Fri Sep 18 13:01:22 EDT --- > Last modified: Wed Mar 3 13:26:47 EST Index:Makefile.ALPHA 9c9 < CFLAGS = -DPVM -D$(PVM_ARCH) -DLITTLE_ENDIAN -I$(FMRI)/include/$(ARCH) -I$(PVM_ROOT)/include -O --- > CFLAGS = -DPVM -D$(PVM_ARCH) -DLITTLE_ENDIAN -DFORTRAN_ADD_UNDERSCORE -I$(FMRI)/include/$(ARCH) -I$(PVM_ROOT)/include -O Index:Makefile.SUN4SOL2 9c9 < CFLAGS = -DPVM -D$(PVM_ARCH) -I$(FMRI)/include -I$(PVM_ROOT)/include -O --- > CFLAGS = -DPVM -D$(PVM_ARCH) -DFORTRAN_ADD_UNDERSCORE -I$(FMRI)/include/$(ARCH) -I$(PVM_ROOT)/include -O Index:Makefile.common 29c29 < partialk $(ARCH_SRCDIRS) --- > partialk physio ireg3d $(ARCH_SRCDIRS) Index:src/brain/Makefile 12,13c12,15 < PKG_LIBS = $(LAPACK_LIBS) -lcrg -lmri -lpar -lbio -lacct -lmisc -larray -lm < PKG_MAKEBINS = $(CB)/vmpfx $(CB)/generate_vmpfx_in $(CB)/byte_extract --- > PKG_LIBS = $(LAPACK_LIBS) -lfmri -lcrg -lmri -lpar -lbio -lacct -lmisc \ > -larray -lm > PKG_MAKEBINS = $(CB)/vmpfx $(CB)/generate_vmpfx_in $(CB)/byte_extract \ > $(CB)/smooth_missing 18c20,21 < drift.c lbfgs.c qnewton-lbfgs.c tmperror.c byte_extract.c --- > drift.c lbfgs.c qnewton-lbfgs.c tmperror.c byte_extract.c \ > smooth_missing.c 21,22c24,25 < README < MISCFILES= vmpfx_proto.t ramps.txt --- > smooth_missing_help.help README > MISCFILES= vmpfx_proto.t vmpfx_null_proto.t ramps.txt 85a89,97 > $(SINGLE_HELP_LD) > > $O/smooth_missing.o: smooth_missing.c > $(CC_RULE) > > $O/smooth_missing_help.o: smooth_missing_help.help > $(HELP_RULE) > > $(CB)/smooth_missing: $O/smooth_missing.o $O/smooth_missing_help.o ${LIBFILES} Index:src/fmri/Makefile 11c11 < PKG_EXPORTS = fmri.h glm.h smoother.h parsesplit.h --- > PKG_EXPORTS = fmri.h glm.h smoother.h parsesplit.h quaternion.h fshrot3d.h 13c13 < PKG_MAKEBINS = $(CB)/smoother_tester --- > PKG_MAKEBINS = $(CB)/smoother_tester $(CB)/quat_tester 18,19c18,20 < smoother_tester.c parsesplit.c fft3d.c < HFILES= fmri.h lapack.h glm.h smoother.h parsesplit.h --- > smoother_tester.c parsesplit.c fft3d.c quaternion.c quat_tester.c > HFILES= fmri.h lapack.h glm.h smoother.h parsesplit.h quaternion.h \ > fshrot3d.h 25c26,27 < $O/smoother.o $O/parsesplit.o $O/fft3d.o --- > $O/smoother.o $O/parsesplit.o $O/fft3d.o $O/quaternion.o \ > $O/fshrot3d.o 28c30 < $O/parsesplit.o $O/fft3d.o --- > $O/parsesplit.o $O/fft3d.o $O/quaternion.o $O/fshrot3d.o 41a44,46 > $O/fshrot3d.o: fshrot3d.c > $(CC_RULE) > 57a63,74 > $(CC_RULE) > > $O/quaternion.o: quaternion.c > $(CC_RULE) > > $O/quat_tester.o: quat_tester.c > $(CC_RULE) > > $(CB)/quat_tester: $O/quat_tester.o $L/libfmri.a > $(SINGLE_LD) > > $O/chirprot.o: chirprot.c Index:src/ireg/Makefile 11c11 < PKG_MAKEBINS = $(CB)/ireg --- > PKG_MAKEBINS = $(CB)/ireg 32a33,36 > > $(CB)/ireg_chirp: $O/ireg.o $O/ireg_help.o $O/../fmri/chirprot.o $(LIBFILES) > $(LD) $(LFLAGS) -o $B/$(@F) $O/ireg.o $O/ireg_help.o \ > $O/../fmri/chirprot.o $(LIBS) Index:src/mri_util/Makefile 15c15,16 < $(CB)/mri_change_dims $(CB)/mri_counter $(CB)/mri_type_convert --- > $(CB)/mri_change_dims $(CB)/mri_counter $(CB)/mri_type_convert \ > $(CB)/mri_pad 24c25,26 < mri_change_dims.c mri_counter.c mri_type_convert.c --- > mri_change_dims.c mri_counter.c mri_type_convert.c \ > mri_pad.c 31c33,34 < mri_counter_help.help mri_type_convert_help.help --- > mri_counter_help.help mri_type_convert_help.help \ > mri_pad_help.help 114a118,126 > $(SINGLE_HELP_LD) > > $O/mri_pad.o: mri_pad.c > $(CC_RULE) > > $O/mri_pad_help.o: mri_pad_help.help > $(HELP_RULE) > > $(CB)/mri_pad: $O/mri_pad.o $O/mri_pad_help.o $(LIBFILES) Index:src/csh/ts_brain.steps.csh 18c18 < echo '# $Id: ts_brain.steps.csh,v 1.6 1998/08/28 20:07:55 welling Exp $' --- > echo '# $Id: ts_brain.steps.csh,v 1.7 1999/03/24 23:46:42 welling Exp $' 36,37c36,38 < # Repeat deghosting step with merged images < ts.deghost2.csh data/$F_BASELINE2_OUTPUT data/$F_DEGHOST2_OUTPUT --- > # Repeat deghosting step with merged images if desired (don't forget > # to change the first parameter of meanc.csh to data/$F_DEGHOST2_OUTPUT > # ts.deghost2.csh data/$F_BASELINE2_OUTPUT data/$F_DEGHOST2_OUTPUT *** Makefile.LINUX Wed Mar 31 14:51:05 1999 --- /SLOTH-U2/welling/fiasco/development/Fiasco4.2/Makefile.LINUX/ Wed Jan 6 14:01:10 1999 *************** *** 0 **** --- 1,15 ---- + # + # LINUX-specific Makefile for FMRI software + # + # + + CC = cc + CFLAGS = -DPVM -D$(PVM_ARCH) -D_BSD_SOURCE -D_XOPEN_SOURCE -DFORTRAN_ADD_UNDERSCORE -I$(FMRI)/include/$(ARCH) -I$(PVM_ROOT)/include + LD = cc + LFLAGS = -L$(PVM_ROOT)/lib/$(PVM_ARCH) + LIBS = -lpvm3 -lgpvm3 + AR = ar + ARFLAGS= -r + RANLIB = echo ranlib not needed on $(ARCH) for + + include Makefile.common Index:src/fmri/polyt.c 21c21 < static char rcsid[] = "$Id: polyt.c,v 1.3 1998/12/03 00:24:03 welling Exp $"; --- > static char rcsid[] = "$Id: polyt.c,v 1.4 1999/04/01 22:03:38 welling Exp $"; 60c60 < static float npar[3]; --- > static float npar[21]; 455c455 < static float npar[3]; --- > static float npar[21];