Implement BAT calculation

[DrSAR]

I think @firasm is working on that already. Questions:

• how precisely do we need to reproduce Kelly's matlab code? (can we get some test data please with references here)
• can we speed up the calculation (should be vectorizable)
• is it worth considering a totally new algorithm (see here for instance: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3097064/ with the coal miner disaster data being of particular interest)

[firasm]

RE: how's development going.

Good - sorry, went to bed early last night.

I'm currently manually testing to see how it works in different cases. Picking certain pixels in the tumour to make sure that when the BAT is a crazy value, it gets set to 0. Also, I'm testing to make sure that when the BAT cannot be identified (no enhancement in pixel), then the BAT Is set to 0.

Here are a few examples, all of HPG in the same tumour

scn = sarpy.Scan('HPGS4H02.6R1/28')

1. Good example of working as it should

2. Working, but probably shouldn't because the SNR is so bad

3. No enhancement, not working and shouldn't be working

4. Weak SNR but works

5. Low SNR but maybe this should work?

[DrSAR]

Don't we need this to be inf rather than 0?

firasm wrote:

RE: how's development going.

Good - sorry, went to bed early last night.

I'm currently manually testing to see how it works in different cases. Picking certain pixels in the tumour to make sure that when the BAT is a crazy value, it gets set to 0. Also, I'm testing to make sure that when the BAT cannot be identified (no enhancement in pixel), then the BAT Is set to 0.

Here are a few examples 5-good https://cloud.githubusercontent.com/assets/2507459/5319045/66de88c0-7c5a-11e4-9404-00e4ab4edf16.png 4-workingbutprobablyshouldnt https://cloud.githubusercontent.com/assets/2507459/5319049/67031b54-7c5a-11e4-8d17-c03e62c2f30d.png 3-no enhancement https://cloud.githubusercontent.com/assets/2507459/5319048/6701612e-7c5a-11e4-9301-3285bd68dcdb.png 2-weakbutworking https://cloud.githubusercontent.com/assets/2507459/5319047/67003966-7c5a-11e4-8df1-690ee0afd4c9.png 1-low snr https://cloud.githubusercontent.com/assets/2507459/5319046/66fced06-7c5a-11e4-99dd-2df3ba6420e0.png

— Reply to this email directly or view it on GitHub https://github.com/DrSAR/SARlabpy/issues/291#issuecomment-65816312.

[firasm]

remember, we decided 0 was better because when you sum 0s, you can ignore them, but infs will have to be accounted for in a special way. Recall the code statement:

cond = where(runningSum(cond1 + cond2 + cond3 + cond4,3) >=3)[0][0]

Anyway, it's easy enough to tweak later if you prefer infs or nans

[DrSAR]

I remember, but in the final BAT it probably makes sense to have it be Inf. A BAT map is not a map of cond but rather a map of indices that indicate when it went from 0 to a finite value. That jump occurs never, hence index=inf.

[firasm]

okay, adjusted.

BAT = numpy.where(cond==0,Inf,cond)
return BAT

I might have to reshape that if it's an n-d array

[firasm]

okay i need some help in vectorizing the convolve function ...

http://docs.scipy.org/doc/numpy/reference/generated/numpy.convolve.html

Zoom?

[firasm]

First BAT map of NecS3 data...

As @DrSAR would say, 'not bloody bad'..

[firasm]

Same map on a slightly more sensible scale (note that currently the colours are just indices and aren't actually in seconds)

[DrSAR]

can you try on some more data: with HPG and with the recent ones. Also - can we have a comparison to the formerly calculated ones (by Kelly). S

firasm wrote:

First BAT map of NecS3 data...

As @DrSAR https://github.com/DrSAR would say, 'not bloody bad'..

necs3hs05 ik1 11_slice3 https://cloud.githubusercontent.com/assets/2507459/5325389/ade2af0a-7c9e-11e4-9e8f-8549b00df074.png

— Reply to this email directly or view it on GitHub https://github.com/DrSAR/SARlabpy/issues/291#issuecomment-65877435.

[firasm]

HPG animal

[DrSAR]

As far as I'm concerned these look great. One could try a median filter (since there is a bit of salt-and-pepper artefact visible). Other than that - can we try histograms on this now? And also a direct comparison with previously obtained maps for known animals?

[firasm]

I was working on comparing this to previous BAT maps and ran into a couple of issues (#289 and #294)

Working on resolving them now..

as for histograms,

I remember you coded a nice adata histogrammer - could you email me a notebook of it in use so I can refresh my memory on how to use it?

[firasm]

Here's the scoop

• BAT is now working on 4D data as well as 3D data (#293)
• BAT now correctly uses the right amount of DCE data to estimate the injection point (#294)
• I've loaded up some of KM's BAT data and am comparing our maps and we're pretty different

see below:

Some issues outstanding that are preventing direct histogram comparison:

1) the code to read bruker data in matlab doesn't properly account for an offset, so the tumour isn't quite at the same place in both images.

2) Since I'm looking at old data, we don't have the ROIs loaded in so I look at the entire image BAT map. Problem here is that because of problem 1, it's not straightforward to just load the ROIs and multiply. to address this, I'll have to figure out the offset from the method file and shift it first.

Doable, but annoying...

3) The overall shape of the tumour seems to be matching (where BAT was correctly calculated), but the BAT in seconds values appear to be offset.

4) This could be due to differences in how we converted indices to times, because I'm doing it on signal intensity (which is inherently less noisy and less prone to fit errors) and KM calculates it on concentration. This code technically works on concentration as well, but because its old data, we don't have LL T1 maps, just the VFA method (and code for calculating T1 using vfa hasn't been written yet) so I can't compare apples to apples.

Need to take a break for a bit...but i'll be back at this later tonight - advice on what to do next would be good. My thought is to try and investigate why our BAT values are offset (by ~ 12 seconds or so) as described in problem 3 above

[firasm]

P.S. The BAT algorithm completely failed on 'HPGS13testMouse.sV1/8' (the new mouse we did last week) because there was no enhancement detected according to the rules we set out.

[DrSAR]

Does the MATLAB code do any shifting or does it simply ignore the issue?

[firasm]

Near as I can tell, this is how the data is read in:

    if ~isempty(dcedirs)
        display(sprintf('DCE data loaded from %2.0d',dcedirs))
%         dce_data = load_bruker_data_and_header(num2str(dcedirs));
        data = read_rare_inv_prep(num2str(dcedirs));
        dce_data.hdr = loadBrukerHeader(num2str(dcedirs));
        for k = 1:size(data,4)

%% THIS LINE
        dce_data.img(:,:,1,k) = circshift(fftn(data(:,:,1,k)), [floor(size(data,1)/2)-1 floor(size(data,2)/2)-1 0]);
%% END LOOKING

        end
        dce_TR = dce_data.hdr.TR;
        dce_data.hdr.time_res = dce_data.hdr.TR*dce_data.hdr.image_dimension(2)*dce_data.hdr.image_dimension(3)/1000;
    end

[DrSAR]

This is not really a shift. It's assuming 0,0,0 offsets. We could apply our shifts to the data and they should line up. There is this issue that our dimensions also have to be aligned. Are they?

firasm wrote:

Near as I can tell, this is how the data is read in:

if ~isempty(dcedirs) display(sprintf('DCE data loaded from %2.0d',dcedirs)) % dce_data = load_bruker_data_and_header(num2str(dcedirs)); data = read_rare_inv_prep(num2str(dcedirs)); dce_data.hdr = loadBrukerHeader(num2str(dcedirs)); for k = 1:size(data,4) dce_data.img(:,:,1,k) = circshift(fftn(data(:,:,1,k)), [floor(size(data,1)/2)-1 floor(size(data,2)/2)-1 0]); end dce_TR = dce_data.hdr.TR; dce_data.hdr.time_res = dce_data.hdr.TR_dce_data.hdr.image_dimension(2)_dce_data.hdr.image_dimension(3)/1000; end

— Reply to this email directly or view it on GitHub https://github.com/DrSAR/SARlabpy/issues/291#issuecomment-66239003.

[firasm]

no, dimensions not lined up - i have to rotate my image CCW 90 degrees to match the BAT Input

[DrSAR]

ours should be aligned with BRUKER images, so you should probably rotate Kelly's images CW 90?

firasm wrote:

no, dimensions not lined up - i have to rotate my image CCW 90 degrees to match the BAT Input

— Reply to this email directly or view it on GitHub https://github.com/DrSAR/SARlabpy/issues/291#issuecomment-66258660.

[firasm]

bumping this up, now that this is again a priority before next week's imaging experiment.

I have created a barebones example in an ipython notebook for you to work on with my attempt at implementing the new BAT collection, and comparing it with Kelly's version (we have now sorted out the dimensions, rotations, etc...)

[firasm]

P.S.this notebook can be found in HPG/new BAT Algorithm

[firasm]

Hey Stefan,

Just a reminder that we need this ready for Friday's HPG experiment. Preferably sooner so we can play with it a bit.

Firas

[firasm]

Current comparison of FM BAT calculation and KM BAT calculation:

Not bad, but I feel like JB would be be upset because my method nicely homogenizes everything and removes the long BATs

[DrSAR]

dare we try a correlation plot or is the pixel-by-pixel alignment not good enough for that?

[firasm]

Alignment should be perfect ish..., so that should be doable

On Wed, Apr 8, 2015 at 8:58 PM, DrSAR notifications@github.com wrote:

dare we try a correlation plot or is the pixel-by-pixel alignment not good enough for that?

S

Reply to this email directly or view it on GitHub: https://github.com/DrSAR/SARlabpy/issues/291#issuecomment-91109641

[DrSAR]

This doesn't look like a 1:1 match (it's also heavily zoomed on the 'reasonable' range of BAT values).

[firasm]

So... This means we're in trouble right ?

The BAT algorithms aren't related by some linear offset, so the algorithms are actually measuring different things.

Other than the exceptions KM coded in, there is one other major difference in the algorithm - KM manually entered in the injection point for each scan whereas I used the full mouse average and an automatic function to estimate the injection point.

Any suggestions on what to do next here ? We'll need this in the afternoon today as we inject HPG.

P.S. Can you be there this afternoon? Around 1:30 pm?