July31T-con/Synopsis

Each of the algorithms has an optimal voxel dimension and they are not the same. We must decide what format of the data will be used as the common starting point.

I have three goals here:

1. to summarize the new information gained during the past 9 days since we had our T-con
   
2. to capture the essence of the discussions that took place in pursuit of that information so that no one is left out inadvertently from the discussion and
   
3. to propose a final solution regarding the format of the starting data to be processed by all the groups for the Santa Fe meeting that will be a reasonable compromise between the in principle optimal but completely impossible solution of all groups using the "same data" and the "every group for themselves" solution that will require more work in Sante Fe to sort out our results.

3. Proposed Solution
To cut to the chase: We will post 2 data sets. The "raw" data is uploaded now. Marc Niethammer will prepare and upload isotropic down sampled data, ROIs and computed tensors. He will provide a clear detailed description of what exactly he does as part of the upload, but see email trail below for his proposal. Those groups that want to work with upsampled isotropic data will generate their own upsampled data sets from the raw data and will include that step in their slide presentation as part of their processing methods. We will discuss in Sante Fe whether we should post the data and the resampled ROIs for the others to use. If you want to also use the eddy current corrected data and then show quantitatively the improvement then do so using the raw data.

1. New Information- Key points extracted from the discussion below

• The wiki has been edited to correct the image acquisition details. The data were collected on a GE scanner using parallel imaging. The prescribed slice resolution is not what is actually written in the "raw" image file as provided by GE. There is some proprietary preprocessing done as or after the data are taken out of K-space and before the non-linear combination of the parallel images. The data available on the scanner has a resolution of 0.9375 x0.9375 x1.7mm (no gap), FOV 240 x 240. That data is the closest to "raw" that is available and is what is currently posted as dwi. The ROIs were generated at this resolution.
• There are three camps within our small group, those whose algorithms can work with this data "as is", those who need upsampled isotropic voxels and those who need downsampled isotropic voxels.
• The dwi-EdCor files are the "raw" data that have been corrected for eddy currents using an FSL tool to do an affine registration of each dwi to the first baseline.
• The ROIs for the Fornix were corrected and replaced on August 3, 2007. Please be sure you are using the correct ones.
• The nrrd headers for the data sets were corrected and replaced on August 2, 2007. Please be sure you are using the correct ones.
• Per Sylvain and Karl, the eddy current corrected images are much better qualitatively. No quantitative metric for this is provided.
• Marek provided the requested information required for clinical relevance of this project WRT the choices of tracts and potential hypotheses to be tested. I have edited them into the ROI page and our Project Page.
• Per Stephan Aylward, Core 2 engineers are interested in helping us to post our final results and/or any other material we wish to share with the scientific community. The top-level website for the effort is http://insight-journal.org/nlm/
• Outstanding questions:

1. What is the b value of the scans?
2. Is the b value included in the header? If not, can it be added?

2. Essence of Email discussion
Attempted to reconstruct several parallel threads, the following are 'mostly' in chronological order. --rlg
"On Jul 31, 2007, at 4:57 PM, Marc Niethammer wrote:

Hi Guido,

it is not completely clear which particular algorithm is implemented on the scanner for the MR signal reconstruction. There is certainly some form of scanner-internal up-sampling, however, to reverse-engineer this may not be trivial.

To make sure that the comparison between the methods is fair, I agree that it is desirable that all of them run with data that is, as far as possible, identical. This will only be possible to a certain extent, since some algorithms work on the diffusion weighted images directly, others use the tensors as their input, but most (all?) of the algorithms perform some form of algorithm-internal data interpolation anyway. (This is easiest to see for the classic streamlining algorithm that needs to compute values with subpixel resolution.)

Here is a possible suggestion of additional datasets we could provide. Let me know what you think.

(1) Downsampled datasets (original and eddy current corrected) to 1.7mm isotropic voxels. The downsampling would be performed on the diffusion weighted images using Gordon's teem library (using a Catmull-Rom interpolation kernel). This includes downsampled labelmaps using nearest neighbor interpolation.

(2) Estimated diffusion tensors (using simple linear estimation) for the datasets in their original resolution as well as for the downsampled datasets. This should be the input data for all methods relying on diffusion tensor input.

Best,

Marc"

On August 1, 2007 Guido wrote: "C-F and all,

thanks for these hints. The problem that I see is that for the upcoming meeting, we might not have a lot of time to carefully evaluate how exactly GE does the upinterpolation. I agree that they are "not stupid" and surely do a smart trick so that the data looks better. The question would be if we really would know.

Marek wrote me that the original voxel resolution is "240 FOV and 144 x 144 in plane resolution, 1.7mm slice thickness and no gap". This would be a native voxel resolution of 1.67mm, but we finally get something like 0.9375mm for a 256 matrix. So it seems that the original matrix of 144x144 is up-interpolated to 256x256 (factor 1.7777.) due to some algorithm, it might be zero-padding or something much more complex, but definitely not factor of 2. Also, the nrrd header says 1.7mm for spacing in z and 2mm for the thickness in z, logically that would imply that we have a gap of 0.3mm in z-direction. However, the z-spacing that our algorithms finally use (reading the nrrd header with ITK readers) for the processing is indeed 1.7mm. Could it be that the z-thickness in the nrrd header is wrong? It would imply an overlap of the 2mm slices with 1.7mm spacing which is a bit unusual. Our nrrd specialists might be able to solve this mystery.

Sorry about all these numbers, but the whole "DTI comparison group" would need to use the same data. We might be more wise after very careful exploration, but we need to make a decision very quickly. Some algorithms have problems with full size data matrix of 256x256x81 voxels with nonisotropic voxel size of 9.9375x0.9375x1.7. Up-interpolation as John at GA does is another solution for getting isotropic voxels, but do we really gain information - such huge data with 256x256x160 matrix also creates computational problems for some of our groups.

I would somehow be inclined, as also C-F suggests, to leave the z-spacing at 1.7mm and downsample the slices back to something like 1.7mm (as they seem to have come from a native measurement of 1.6666mm). The question to all of the signal processing experts is if we could just use SINC interpolation for the 2D slices without "damaging" the information too much. I think we should decide relatively soon to get the datasets ready for all the groups.

Btw, there are the raw data and also the EPI corrected data (using FSL) on the NAMIC download page. It seems there is so far not yet a quantitative analysis of how much the EPI correction improves the quality. Due to our time pressure, I would also suggest that we could put this issue on hold and just use the uncorrected data for our first 10 case comparison study.

This issue of native resolution versus scanner interpolation etc. could be a good topic to be added to C-Fs recommendation to talk about linear/nonlinear pre-filtering and interpolation issues since these change the data significantly.

Guido

Carl-Fredrik Westin wrote: Regarding voxel size I agree that we should try to get back down to the original voxelsize. If the image data upsampled with zero padding in the Fourier domain, which can be checked by taking the Fourier transform and see of the signal is zero on the upper half of the spectrum. If so the zeros can be removed by cropping the Fourier data be and transformed back to the original size.

Unfortunately I think there is more to it, sinc interpolation that zero padding implies can give nasty ringing artifacts, and some additional filtering can reduce that, and GE has probably figured that out too. In any case, we should try and look at the result. In the end I think we will do little damage to the data by down sampling it back to the original size, independent of the method. Simplest method is to pick every second sample, but only works if the up-sampling factor is two.

A side note, for the Riemannian interpolation community we should probably remove any linear pre-interpolation, which zero padding is.

C-F

-- Dr. Carl-Fredrik Westin Director, Laboratory of Mathematics in Imaging (LMI) Associate Professor of Radiology, Harvard Medical School Brigham and Women's Hospital, 75 Francis St., Boston, MA 02115

Phone: (+1) 617-525-6209, Fax: (+1) 617-582-6033 Email: westin@bwh.harvard.edu, Office: 1249 Boylston St, 2nd floor http://lmi.bwh.harvard.edu/~westin

On Jul 31, 2007, at 7:26 PM, Marek kubicki wrote:

Hi Guido, Sorry I could not participate in the t-con today, but we had 5 subjects scheduled for DTI and fMRI scans today, and I had to run the scanner. I remember having the same discussion with CF and our GE engineer some time ago. Originally data is acquired with 240 FOV and 144 x 144 in plane resolution, 1.7mm slice thickness and no gap. We do not preprocess the data (besides eddy current distortion correction, which makes them less blurry), so the original data is by default upsampled to higher resolution on the scanner, by GE software, and our engeneer did not know if there was the way to avoid this. We were thinking that it makes sense to downsample it to the original resolution, which was supposed to be 1.7 mm cube... Other thoughts? Marek"

"On Aug 2, 2007, at 1:27 PM, Marek Kubicki wrote:

Hi Guys, I looked at the action items from the t-con, and I have some thoughts/comments/questions for all of you, regarding my part: 1. I edited some information regarding the data- FOV is 240 mm square, not 140. We have eddy current corrected, not distortion corrected images. We do not have distortion correction tools. For the eddy current distortion correction, we took all the gradient images, and registered them together. I think its linear registration, nothing fancy. 2. I will try to get the information from GE, about how they up-sample the images, and if there is anything that can be done to reverse the process. 3. Regarding the clinical relevance, this is tough. We picked tracts that connect frontal and temporal lobes, and all of them have been shown, at some point, to be abnormal (FA lower) in schizophrenia. I will put a slide that illustrates their exact function, and why we think they are affected in schizophrenia, if this is what you are looking for, but there is not much that has been done in terms of specificity of the white matter abnormalities. So it might be tough to directly answer to your question- what we expect to find in schizophrenia on this dataset. We have very small datasets in this project, so I would not expect to find statistical group differences. Because of the huge symptom variability within the schizophrenia spectrum, we usually need at least 20 cases, plus the same number of controls, to show the difference. I guess we could look at the effect sizes, and see how many subjects would we want for each method and each tract to show group differences? We could also look at the measurement variability within the control group? Or how combination of measurements increases the sensitivity of the DTI measurements? Your thoughts? Marek"

"Randy and all,

I followed the email discussion but did not yet step in. Good that we solved the 1.7mm slice thickness, this is what we all thought and used, and those who often see brain images would notice that using 2mm, the brains would look unusually scaled in z-direction.

Casey did a Fourier transform with Matlab (see attachment) and we conclude that there is a real cutoff of the signal at the frequencies which originate from the 144x144 grid sized image. It is clear that this 144x144 got interpolated to 2556x2556 by GE on the scanner. What we propose is that there is no harm to downinterpolate the 256x256 back to 144x144 (or voxel size of 0.9375mm2 to 1.67777 (240FOV/144)). Casey proposes a "good" interpolation, like a sinc approximating kernel 9windowed since or high-order b-spline).on the 2D slices (2D only!) to get to a matrix of 144x144x81, these interpolations are available in the ITK toolkit. WHO WOULD VOLUNTEER TO DO IT?

Second, I agree that for the coming meeting, we might use the downinterpolated, non-corrected data (and as Randy says everybody should feel free to spend the effort to quantitatively show the improvement of the EPI correction, which I believe should be measurable would we know how exactly to compare pre/post-correction). As I suggested earlier, the precorrection is a discussion point along the lines of other preprocessing discussion proposed by C-F, like prefiltering the data with different choice of intepolation schemes.

Guido

Randy Gollub wrote: Ahhhh, DICOM header clarity! Yes, please remove the incorrect nrrd-headers and replace with correct ones. This resolves one key aspect of the discussion. Now what about up and down- sampling? And are we all in agreement that everyone should work from the non- eddy current corrected data set to start with?

Perhaps those who have time can also use the eddy current corrected version and then present their results of the difference it made. It could be an incremental step forward to boot.

Randy

On Aug 2, 2007, at 2:27 PM, Marc Niethammer wrote:

Hi everyone!

According to the DICOM header (see below) the thickness is indeed 1.7mm. So there is no gap. I can upload updated nrrd-headers if desired.

Marc

pnl_d6_4:/projects/schiz/diff/parallel_diffusion/caseD00959/000004.SER% /projects/mark/bin/print_header 000001.IMA print_header V1.0 cross_product = (-0.001808,0.031548,0.999500) cross_product(abs) = (0.000000,0.000000,-1844237796292823455052365436477544292218359422354176133202280781869652014465647948127420522205374599172130974934072049841985697262590661287504171604736610899861487636112302173053691298350732573252851412225942032239655277432678131009512974993123543133995413985171543706832834410864132103366101836447403212800.000000) axial! input_prefix = 000001 output_prefix = suffix = <none> filename_format = %s.%03d patient_name = patient_sex = date = 20061107 time_hr_min_sec = study_desc = MR085 MRI HEAD (SCH) series_desc = TENSOR ASSET hospital_name = Brigham and Womens 3 Tesla patient_id = exam_modality = MR image_type_text = dicom image_type = ORIGINAL\PRIMARY\OTHER sop_instance_uid = 1.2.840.113619.2.176.3596.6864640.8898.1162906008.899 image_time = 150444 manufacturer = GE MEDICAL SYSTEMS referring_physician = physician_reading_study = N/A manufacturers_model_name = SIGNA EXCITE anatomic_structure = N/A birth_date = weight = body_part_examined = N/A repetition_time = 17000 echo_time = 74.1 inversion_time = 0 A reconstruction_diameter = 240 generator_power = N/A series_in_study = N/A images_in_acquisition = 4602 rescale_intercept = N/A rescale_slope = N/A requested_proc_desc = N/A patient_age = 35 x_resolution = 256 y_resolution = 256 bytes_per_slice = 131072 header_size = 10846 slice_number = 0 image_type_num = 6 byte_order = 1 status = 0 number_echoes = 0 echo_number = 1 compressed = 0 first_slice = 0 last_slice = 0 num_missing = 0 bytes_per_pixel = 2 num_bytes_data = 0 num_bytes_header = 0 bits_stored = 60904 high_bit = 60904 bits_allocated = 60904 patient_position = HFS patient_orientation = N/A number_of_slices = 0 exam_number = 4102 series_number = 4 gantry_tilt = 0.000000 pixel_xsize = 0.937500 pixel_ysize = 0.937500 fov = 240.000000 aspect = 1.813333 thick = 1.700000 space = 0.000000 image_location = -38.490765 coord_center_r = 5.762001 coord_center_a = 26.738998 coord_center_s = 84.967499 coord_normal_r = 0.000000 coord_normal_a = 1.000000 coord_normal_s = 0.000000 coord_r_top_left = 125.762001 coord_a_top_left = 146.738998 coord_s_top_left = -35.032501 coord_r_top_right = -114.237999 coord_a_top_right = 146.738998 coord_s_top_right = -35.032501 coord_r_bottom_right = -114.237999 coord_a_bottom_right = -93.261002 coord_s_bottom_right = -35.032501 imagepositionpatient = -125.762001 -146.738998 -35.032501 imageorientationpatient = 0.999710 -0.023926 0.002564 0.023988 0.999216 -0.031495 pnl_d6_4:/projects/schiz/diff/parallel_diffusion/caseD00959/000004.SER%

Hi Everyone, This may have been resolved already on traffic that I haven't seen, but was there, in fact, a slice gap in this data? That seems a little odd for DTI data that you're going to use to generate tracks. You could probably figure out whether or not this was the case from the absolute slice positions that are sometimes stored in the DICOM header.

Karl"

On August 3 Marek wrote: "More info from GE:

1) Zero filling is applied first, followed by the FFT, followed by the parallel imaging processing - SENSE

2) A Fermi filter is applied to the raw data which rolls off the higher spatial frequences. The Fermi filter can be turned off with a cv. The Fermi filter is applied to attenuate ringing which MIGHT occure due to the finite sampling of the acquired data. The Fermi filter has a value of 1.0 at the center of k-space and roles off expenentially to 0.5 at a radius of yres/2 and xres/2. If the noise dominates the k-space data at the edges, then the Fermi filter has little or no effect. This is a function of RF coil, SRN, scan parameters (slice thick, FOV, yres, xres, etc). The fermi filter has a window and radius which determins how fast the filter roles off.

3) Zero filling does not add artifact to the images. In simply provides higher resolution to what exists. There is no information within Zeros. Any edge or ringing will be more apparent. It also does not affect SNR.

---

Original Message -----

From: Carl-Fredrik Westin To: Marc Niethammer Cc: Karl Helmer ; Randy Gollub ; Guido Gerig SCI Utah ; Marek Kubicki ; Sylvain Bouix ; Guido Gerig ; fletcher@sci.utah.edu Sent: Friday, August 03, 2007 1:24 PM Subject: Re: REMINDER: NA-MIC Tractography Measures Conference T-CON tomorrow

We just wanted to figure out if we could undo what GE did. Since there are some, although weak, energy in the higher part of the spectrum, something non-linear must have happened in addition to the zero filling. Only an addition linear filter would not invent any new frequencies

Does anyone know if the parallel reconstruction is happening before or after the zero filling? If it the signal is only from quantization, I think we would see some edge to the higher frequency area. It could well be that the image are zero-filled first, transformed to the spatial domain, and then reconstructed using the parallel machinery.

In any case, I think the conclusion is that we should not force any specific down sampling scheme for everyone. Whoever wants to do it can make it a part of their algorithm. C-F"
On Aug 3, 2007, at 1:10 PM, Marc Niethammer wrote:

Hi Karl!

No, I haven't tried this, since I won't be able to compare results so easily then. I just wanted to see if I can go back and forth between spatial and Fourier domains. If I could, we could do whatever we would want based on the Fourier data. It seems like we cannot though (or at least not in a very straightforward manner). If there is some form of scanner discretization (Is there?) we won't be able to perfectly reconstruct. And my feeling is that there must be, because for a Fourier-based reconstruction the results will certainly be real, whereas the data delivered by the scanner is integer data.

Marc

Hi Marc, So you essentially re-zero-filled already zero-filled data and then did the IFT back into the zero-filled matrix size. Did you try putting the LP-filtered data into a 144 x 144 matrix and doing the IFT to see what that looks like? That difference I'm betting will be bigger.

Karl

Hi everybody,

I played a little bit with the diffusion data. The results are here

http://www.na-mic.org/Wiki/index.php/July31T-con/PlayingWithTheDiffusionData

Marc"

From Guido on August 6: "Randy and all,

I did not update the Wiki, and this is indeed a very interesting discussion. Although C-F proposes that we should not force anyone to accept a specific data interpolation, I fear that any comparison a the forthcoming meeting becomes very difficult. I understand the concerns of C-F and Marc that cutting the high frequencies indeed seems to cut some information. Would some labs us upinterpolated, no-interpolated, down-interpolated data with different interpolation schemes, everyone would work with different input data which by nature already would be seen as systematic differences. The telephone conference clearly showed that some groups really can't use the raw data but would down-sample them to an isotropic grid. Forcing all the participating groups to use this non-isotropic data might "cut" some labs out of this comparison. Maybe we really have to make a compromise just for this comparison to provide a downinterpolated, close-to-isotropic dataset for everyone and leave the more sophisticated analysis of optimal "redoing of the GE upsampling" or the use of the raw GE data for additional research to be discussed at the meeting?

Randy and all, could we propose the following: a) we decide about using the non-EPI corrected data for the time being (given the short time left for the meeting, not because EPI correction might not be useful but because we first need to quantitatively show its advantage)

b) we provide downsampled images (sampling the 256x256 back to a 144x144 grid with the standard procedure as used by Marc and provide this close to isotropic dataset to everyone. This will ensure that every group has access to a standard, close to isotropic dataset (1,67x1.67x1.7).

c) leave it to every group to use the original matrix, more sophisticated down- or upsampling, and to compare - and discuss implications/differences as part of the report.

Please understand that I don't want to stop the current, exciting discussion about the optimal reformatting of the data or redoing the GE upsampling, but I would like to see this as a component of the workshop itself. I think that all these questions is a very important part of the NAMIC DTI toolkit and recommendations for users, and of course providing the tools to perform a correct/optimal data preprocessing."

Later that day Sylvain expressed the following concerns: 1) All ROIs were drawn on the original data and they will have to be downsampled too. 2) "WRT to Guido's point > c) leave it to every group to use the original matrix, more > sophisticated down- or upsampling, and to compare - and discuss > implications/differences as part of the report. I am somewhat less convinced about this. I say if we want to be consistent, then only one single dataset should be provided so that all methods start with the same data.

Unless there is *strong* opposition against this, the PNL will prepare for each case: 1- A downsampled 1.67x1.67x1.7mm data set of DWIs in NRRD format. 2- Its corresponding plain vanilla linear least square tensor fit in NRRD format. 3- The associated ROIs downsampled to same resolution in NRRD format. All other data sets should be removed from this project (original and epi corrected).

For consistency, if the candidate tracking methods uses DTI as input then the tensor data *must* be used. Differences in tensor estimation techniques should not justify poor tracking, or at least should be observable in all the tracking techniques. Similarly, If the candidate tracking methods works directly on the DWIs then the DWI data *must* be used.

There is no perfect solution to finding the ideal data set. This is the lowest common denominator to all techniques within the NAMIC community and beyond. Given the large variability in processing techniques the least we should do is use the exact same input.

Best, -Sylvain

Return to July31T-con Page