ARTREPAIR INSTRUCTIONS --------------------------------------------------------------- This protocol provides a visual review and runs semi-automatic detection and repair functions. These programs do not erase files, so no harm can be done to your raw data. The programs will write new files to the same directory if repairs were required. NOTE: Sample figures and outputs are in the document ArtRepairOverview.ppt. CONTENTS 5.1 CONTRAST MOVIE to review data 5.2 NOISE FILTERING to remove and reduce slice and spike noise 5.3 ARTIFACT REPAIR to repair volume artifacts 5.4 REPAIR AND COMPARE - one step method to repair and reanalyze existing SPM design 5.5 GLOBAL QUALITY to validate estimation accuracy 5.6 ADVANCED FUNCTIONS - Multiple session experiments, sample batch scripts ----------------------------------------------------------------------------- 5.1 CONTRAST MOVIE ----------------------------------------------------------------------------- This program allows a user to view every individual voxel in every volume. For a quick review, one could run art_global (Section 5.3) to see where the most problematic scans are located, and just use the movie on those scans. Start the movie/slider program by typing >> 'art_movie' at the Matlab prompt. Type 'help art_movie' for online help. INPUTS Choose Orientation: Transverse, Sagittal, or Coronal. A * marks the orientation of slice collection which gives the best views. Choose a Set of Images (in .img format) to review. Choose a Range of images to process, e.g. 51:100 will show 50 volumes, starting at 51. The program runs well with 50 or 100 images depending on your computer. Hundreds of images can be chosen, but the program will run slower. Select All slices, or 25 close-up slices. In close-up mode, about 20-25 consecutive slices are shown in a montage at twice the image size per slice. Works well for raw images. In All Slices, every slice is shown. Works well for normalized images. Select Raw or Contrast or High Contrast display. Raw shows the image data in the set of scans. Contrast Mode is the best for artifacts. The display shows yellow as higher, and blue as lower from the reference image. Values are amplified to be 5x more sensitive for viewing. Contrast mode helps to see artifacts in raw images. High Contrast is the same as Contrast, but amplified 20x instead of 5x. Allows a close look at voxel noise. Select Reference image. User may choose any image, or allow the the program default to the second image in the range. Typical user choice might be a mean image or a rest state. Select Movie & Time History, or interactive Slider. Movie mode allows a rapid review of all the data. If a movie, select a movie frame rate. Three loops are played. 1 or 2 fps is slow enough to visually spot artifacts. At the end of the movie, a time history of the average intensity and position of each scan is displayed. Slider mode allows interactive views of individual volumes and inspection of individual voxels using zoom in and out. Select Export option, if desired, for Analyze image output. REVIEW THE DATA The contrast mode uses a colormap scale of 160, where bright yellow or bright blue represent positive or negative deviations from a reference image. The internal scaling works such that the max raw image intensity is about 2000 (dimensionless units). For fMRI data, the brightest signal is usually CSF, and gray matter is dimmer. Thus, if the gray matter signal after scaling is 1600, the contrast colormap range shows 10% deviations. Since activation sizes are typically 2%, a contrast image of good data should just look dark. Any bright colors are probably artifacts. The powerpoint slides have examples of normal and abnormal data. ----------------------------------------------------------------------------- 5.2 DETECT AND REPAIR SLICE ARTIFACTS (NOISE FILTERING) ----------------------------------------------------------------------------- The bad slices show up best in the raw images, and the bad slice data spreads around in slicetiming and realignment. So check for bad slices before any SPM preprocessing. One mode of this program will automatically detect bad slices in the data, repair them with interpolation (from the before and after scans) and write a logfile of the corrections. Start the bad slice detect and repair program by typing 'art_slice' at the Matlab prompt. Type 'help art_slice' for online help. INPUTs Select images. Option: Select Repair Bad Slices and Write Bad Slice Log This program automatically screens the data for bad slices, and will write a logfile showing all the bad slices detected. The slice repaired data is written with a prefix "g". A question pops up about setting a threshold. Leave the default value in place, and adjust it later if necessary according the results of the BadSliceLog (see below.) Option: Median Filter All Data If TR<= 2 sec, you have the option of median filtering the data. This filter is a temporal filter for outliers, and it often helps clean up the data. It is especially good if bad slices were detected, because there are usually moderately bad slices just below the bad_slice detection threshold. The median filter does not write a bad slice log- it just cleans up the data. If you ran the previous option to detect bad slices, we suggest you delete the "g" prefix files from the previous step to save space. The median filter should be run on the original data, and a prefix "f" will be prepended to the name. BAD SLICE LOG If the logfile shows clean data, there were no bad slices detected at the default threshold. Otherwise, the log lists all the bad slices which were corrected, and the percentage of total slices which were corrected. The bad slice repair can hurt activations, so a rule of thumb is to correct fewer than 5% of the slices. If too many slices are corrected, try raising the threshold. If you want to correct a bad slice seen using art_movie, try lowering the threshold. OUTPUT FILES If repairs were done, the output files may have prefixes of "f", "g", and/or "h" depending on the repair method used. ----------------------------------------------------------------------------- 5.3 DETECT AND REPAIR VOLUME ARTIFACTS ----------------------------------------------------------------------------- The volume artifact program works on both global intensity and scan-to-scan movement, so it must be run after realignment. A good opportunity is just before estimation. The input images must be resliced. The program pops up a GUI with suggested scans to be repaired, and suggested scans to be deweighted during SPM Estimation. The GUI is interactive and the user can adjust thresholds. Nothing is written until the user selects REPAIR, at which point a new set of repaired images is written. All the old images are preserved. Start the volume detect and repair program by typing 'art_global' at the Matlab prompt. Type 'help art_global' for online help. INPUTS How many sessions? Usually choose 1. Program allows the first scan of each session to be dropped, if multiple sessions. Select the automatically generated mask The global mean needs a reasonable subject-specific mask. You can build one by hand and use it, but art_automask program works well (at least on the Lucas spiral data.) Have realignment files?. Choose No. For raw data, realignment parameters are not available. The program can show realignment parameters if they've been calculated already. Select image set The program can handle more than a 1000 scans, and from multiple folders. Drop first scan of each session? Sometimes the first volume is a warm-up, or a scanner calibration, and it should not be included in the analysis. Let the program run until the Global Summary Figure pops up. GLOBAL MEAN and SCAN-TO-SCAN MOTION SUMMARY FIGURE The summary figure has four graphs. The top graph shows the time history of the mean intensity within the head region. The time history usually has a slight downward trend, and a good run has only small bumps (less than 1.5% of the mean value). The vertical scale is the dimensionless-units that come off the scanner. The head region is defined as the voxels in the ArtifactMask which was generated during the art_global program. (Use SPM Display function to review the mask ArtifactMask.img.) The second graph shows the size of the bumps relative to the mean of the run. The vertical scale is standard deviations away from the mean. A default threshold line is plotted at 1.5% variation from the mean. The default threshold of 1.5% from the mean is an estimate. For data from the Lucas scanner, this value attempts to keep all data within normal physiological variation and discard outliers due to sudden motion or other artifact. (You can use art_movie to inspect an outlier volume.) A good run has all its peaks below this threshold. The third graph shows the six SPM realignment parameters. (Before realignment is done, this graph shows a crude estimate of the 3 translational estimates of centroid motion. Not to be trusted, it is just locate very peculiar scans.) The fourth graph shows the scan-to-scan movement. The vertical scale is mm/TR. (Rotation assumes a voxel is 80 mm from the origin.) A default threshold line is plotted at 0.3 mm/TR variation. A still subject will have less than this scan-to-scan motion. REPAIR and DEWEIGHT DISPLAY The programs default thresholds are used to estimate outlier volumes and additional volumes to be deweighted during estimation. Outlier volumes to be repaired are specified by red vertical lines. Outliers may be from exceeding the intensity variation threshold, or exceeding the scan-to-scan motion threshold, or by user selection. All repaired volumes are recommended for deweighting. Additional volumes recommended for deweighting are specified by green vertical lines. These volumes will be close to large discontinuities in global intensity or movement EVEN AFTER REPAIR, and so may degrade the accuracy of SPM estimation. Nothing is changed until the REPAIR button is pushed. THRESHOLD ADJUSTMENT BUTTONS The Up and Down buttons raise and lower the thresholds (the changes are linked, no independent changes of each threshold are possible). After each threshold change, the the outliers for repair are recalculated and shown in the figure. OUTLIER EDITING The outlier list can be edited by hand. After Enter is pressed, the outliers for repair are recalculated and shown in the figure. MARGIN BUTTON FOR DEWEIGHTING At any point, hitting the MARGIN button will have the program recalculate additional suggested scans for deweighting. For example, to see recommended deweighted scans, this button could be pressed after any threshold or outlier editing change. REPAIR BUTTON The Repair button will repair any volumes that are outliers at the current threshold. No changes are made until the REPAIR button is hit. The volume repaired data is written with a prefix 'v'. Note that even unchanged images are copied to the new name. The repaired volumes are listed in the text file art_repaired.txt in the same folder. All volumes recommended for deweighting are listed in the text file art_deweighted.txt. (Note the volume index in these lists corresponds to order of the images...the index does not represent the volume number itself unless the first image was V001.) We recommend using the INTERP button as the best method of Repair. This function fills values using linear interpolation from the nearest unrepaired scans. Mean fills values with the mean scan of the run. Despike fills values as a linear interpolation of the immediately preceding and following scans, whether they are repaired or not. OUTPUT FILES If repairs were done, the output files will have prefix "v" added to the name. depending on the repair method used. New files art_repaired.txt and art_deweighted.txt are written to the images folder. The repair list logs which scans were repaired. The deweighted list is for the Repair and Compare program. These files only show the sequence position of the images, so they are NOT accurate if scans are discarded as part of the estimation process. To see the effect of repairs, you can run contrast movie or art_global again on the repaired volumes. ------------------------------------------------------------------------------ 5.4 REPAIR AND COMPARE ------------------------------------------------------------------------------ This is a one step repair and retest process that allows a user to try repairs on an existing result to try to improve the accuracy. All original results and the SPM.mat file are always preserved. Accuracy is measured by a GlobalQuality metric of the range of contrast estimates over the ensemble of all the voxels in the head. The metric assumes that true cognitive contrasts are small relative to artifact and motion noise, so that a reduction in the range of contrasts is evidence that the SPM estimates have become more statistically efficient. Program will copy the existing SPM design to a new folder, [optionally] repair the data with the default thresholds in art_global, replace the input images in the new SPM.mat with the repaired images, and run SPM estimation with deweighting on the repaired images. Global Quality will be measured to check estimation accuracy before and after repair. A comparison summary of accuracies will be shown in the Matlab command window and written to GlobalQuality.txt.(see below). Results depend on the size of the mask image, which is shown in that text file. For the MOST FAIR COMPARISON, run the Global Quality function for each result using the SAME image mask in both runs. If the accuracy after repairs is not good enough (e.g. worse than estimates with the unrepaired data), try a manual repair using the Artifact Repair button (or equivalently >>art_global) and manually adjust the repair thresholds. Then run this >>art_redo program again using the repaired data with the manually adjusted threshold to check the accuracy. Depending on the quality of the data, it will need few or no repairs, or many repairs (even half the data in very bad cases). In some cases, Repairs may make the results worse (in which case, stick to the original result!). A rule of thumb is not to deweight more than 1/3 to 1/2 the data. If the default threshold is discarding too much data, raising the threshold may improve the results. Deweighting is automatically applied in the art_redo function, but it is not directly supported in the SPM manual interface. A sample program art_batchestimation.m shows how the art_deweighted.txt file can be used in batch processing. TROUBLESHOOTING --------------- If one of the Global Quality graphs shows a flat line, it is likely that the mask.img in the original results folder is not a superset of the mask.img in the repaired folder. Use the Global Quality button to run the repaired results with the mask.img in its own folder. (The Repair and Compare button tries to use the same mask for both images to make the fairest comparison.) "Error: Must enter same number of scans as original." A possible cause is if the batch script which generated the original results deleted a few scans (e.g. the first three), but those scans are still in the images folder. Then when art_global wrote v-files for every scan, the number of v-files was larger than the number used in the design. One workaround is to do a manual repair (5.3 Artifact Repair), discard the v-files corresponding to the deleted scans, and then run Repair and Compare using the existing v-files. Another workaround is to rerun the original results without deleting any scans, but take care that the experiment onsets remain OK. ------------------------------------------------------------------------------ 5.5 GLOBAL QUALITY ------------------------------------------------------------------------------ Use the art_summary program to see the effect of repairs on estimation accuracy before and after the repairs are applied. The art_summary program summarizes the global distribution of estimates and ResMS values produced by the SPM estimation process. These summaries are histograms and statistics taken over the ensemble of voxels within the head mask. Generally, for good contrast estimates, the mean of the contrasts should be near zero, the distribution of the contrasts should be small, and the mean of the ResMS should be small. The program draws a figure of the histogram of estimates, and writes a file named "GlobalQuality,txt" to the folder with the input estimates. UNITS ----- The units of Global Quality are PERCENT SIGNAL CHANGE when the peak of the design regressor is 1.0. This choice for the definition allows a user to quickly see the accuracy of estimates relative to the expected effect size of the contrast. Note if the peak value of the design regressor is p (instead of one), the actual percent signal change should be calculated by multiplying the Global Quality measure by p. This peak value can be observed using the SPM Review Design and Explore options. For ResMS, the relative sizes can be used to compare the model residuals, but the units are not percent signal change. INTREPRETING THE GLOBALQUALITY.TXT FILE --------------------------------------- A text file called Global Quality.txt will be written into the Subject folder that includes the results of both GlobalQuality runs. Example Output: GlobalQuality.txt C:\fraX\fr1369\ResultsSPM\con_0010.img Voxels/1000 Mean Std RMS Trimmean 90ile %Vox > 1% AbsMax 13.5820 -0.0090 0.7744 0.7744 -0.0164 0.8316 7.4216 6.1105 9.1010 -0.0024 0.6234 0.6234 -0.0219 0.6789 4.9775 6.1105 4.4810 -0.0225 1.0140 1.0143 -0.0163 1.1520 12.3856 5.8479 C:\fraX\fr1369\ResultsSPM\ResMS.img Voxels/1000 Mean Std RMS Trimmean 90ile %Vox > 1% AbsMax 13.5820 28.5143 15.4679 32.4395 26.9395 22.2822 100.0000 172.9981 9.1010 23.1666 10.4969 25.4338 22.1512 14.0100 100.0000 134.6867 4.4810 39.3756 18.0386 43.3108 38.0512 26.9696 100.0000 172.9981 C:\fraX\fr1369\ResultsRepaired\con_0010.img Voxels/1000 Mean Std RMS Trimmean 90ile %Vox > 1% AbsMax 13.5820 0.0418 0.5663 0.5678 0.0341 0.6163 4.4692 4.0296 9.1010 0.0390 0.4350 0.4368 0.0342 0.4772 2.0108 3.9293 4.4810 0.0475 0.7666 0.7681 0.0360 0.9294 9.4622 4.0296 C:\fraX\fr1369\ResultsRepaired\ResMS.img Voxels/1000 Mean Std RMS Trimmean 90ile %Vox > 1% AbsMax 13.5820 17.3786 10.5526 20.3316 16.1886 15.2464 100.0000 124.2562 9.1010 13.7456 7.0707 15.4576 12.9093 9.0196 100.0000 88.1768 4.4810 24.7574 12.4399 27.7070 23.7724 18.8811 100.0000 124.2562 The key statistic is the Std value in the top row under the con_ images. In the data above, that value is 0.7744 without repairs, and 0.5663 after repairs. Thus, the repairs reduced the STD of the estimation error for this contrast by (0.7744-0.5663)/0.7744 = 27%. Our validation tests show that this value tracks the actual RMS error of test injection examples very well, so the repair was very effective for this case. Another interesting statistic is the Mean value in the top row under the ResMS image. This value is the average (over the head) of the ResMS image from SPM, and thus represents an "average" error in cases when GLM assumptions are valid. In the data below, that value is 28.5143 without repairs, and 17.3786 after repairs. With this measure, the repairs reduced the average ResMS error by (28.5143-17.3786)/28.5143 = 39%. The other statistics that are printed out include: For each image, the first row has statistics for all the voxels, the second row is the inner voxels, and last row is the outer voxels. Errors for the inner voxels are less than for the outer voxels, and results for all voxels are between the inner and outer results. Voxels/1000 is the number of voxels divided by 1000. The statistics are computed over the ensemble of voxels in the head. Values under the con_ will be in percent signal change (when peak of regressor is one.) Values under the ResMS image are relative values, related to tr(RV) from SPM. Mean and Std are the usual statistics. RMS is the combination of mean and Std, as if the truth were zero. Trimmean is the trimmed mean after clipping 5% of outliers on each extreme. 90ile is distance of the 90th percentile value to the trimmean. %Vox > 1% shows the fraction of voxels with extreme contrasts. AbsMax is the most extreme value (min or max) in the distribution. ------------------------------------ 5.6 ADVANCED FUNCTIONS ------------------------------------ Notes for programmers, but these functions have not been tested very thoroughly. MULTIPLE SESSIONS for a subject ------------------------------------ For experiments with multiple sessions, the repair code assumes that all sessions are realigned to the first volume of the first session. The ArtRepair detection process will (or at least is supposed to) repair the scans near session boundaries if the subject moves between sessions. It's assumed that SPM writes the realignment files into each session folder, and that art_global will write the list of files to be repaired into the first session folder. Repaired files are written session by session into separate folders. (art_repair uses the image's file path) Art_repaired.txt for all concatenated sessions is written in the first session folder. (art_repair writes text file to path of first image of the experiment) SAMPLE BATCH SCRIPTS ------------------------------------- Added sample batch script codes for repair and quality check for single session experiments ONLY. art_batchmultisubject Runs multiple subjects through estimation and repair. art_batchrepairsubject For one subject, runs repair, estimation, and global quality. Does art_global, estimation, and art_summary. art_batchestimation For one subject, runs estimation and contrasts with deweighting (Program was named art_batchexample in the earlier version.) art_global Batch call prints jpeg image where title is subject name