%       % Download and unzip energy volume file for gene Rasd2 coronal SectionDataSet 73636089
    mkdir('Rasd2_73636089');
    urlwrite('http://api.brain-map.org/grid_data/download/74819249?include=density', 'temp.zip');
    unzip('temp.zip','Rasd2_73636089');  
    % Download and unzip density volume file for BLAa injection SectionDataSet 113144533
    mkdir('BLAa_113144533');
    urlwrite('http://api.brain-map.org/grid_data/download/113144533?include=density', 'temp.zip');
    unzip('temp.zip','BLAa_113144533');
 
    % Gene expression grids are at 200 micron resolution.
    geneGridSize = [67 41 58];
    fid = fopen('Rasd2_73636089/density.raw', 'r', 'l'  );
    Rasd2 = fread( fid, prod(geneGridSize), 'float' );
    fclose(fid);
    Rasd2 = reshape( Rasd2, geneGridSize );
 
    % Projection grids are at 100 micron resolution
    projectionGridSize = [133 81 115];
    fid = fopen('BLAa_113144533/density.raw', 'r', 'l'  );
    BLAa = fread( fid, prod(projectionGridSize), 'float' );
    fclose(fid);
    BLAa = reshape( BLAa, projectionGridSize );
 
    % Upsample gene expression grid to same dimension as projection grid using linear interpolation
    [xi,yi,zi] = meshgrid(1:0.5:41,1:0.5:67,1:0.5:58); %note: matlab transposes x-y
    d = Rasd2; d(d<0) = 0; % fill in missing data as zeroes
    Rasd2_100 = interp3(d ,xi,yi,zi,'linear');
 
    % Handle "no data" (-1) voxels.
    % Create a mask of "data" vs "no data" voxels and apply linear interpolation
    m = (Rasd2  >= 0); mi = interp3(m,xi,yi,zi,'linear');
 
    % Normalize data by dividing by interpolated mask. Assign value of "-1" to "no data" voxels.
    Rasd2_100 = Rasd2_100 ./ mi;
    Rasd2_100( mi <= 0 ) = -1;
 
    % Create a merged image of one coronal plane;
    gimg = squeeze(Rasd2_100(52,:,:)); gimg = max(0,gimg); gimg = gimg / 0.025; gimg = min(1,gimg);
    pimg = squeeze(BLAa(52,:,:)); pimg = max(0,pimg); pimg = pimg / 0.8; pimg = min(1,pimg);
    rgb = zeros([size(gimg),3]); rgb(:,:,1) = gimg; rgb(:,:,2) = pimg;
    figure; image(rgb);