Working with Aperio SVS files in Matlab – Converting Annotations to Binary Masks

One of the main purposes of having a digital format is to allow experts (e.g., pathologists) to annotate certain structures in the images. Be it nuclei, epithelium/stroma regions, tumor/non-tumor tissue etc. This is easily done with ImageScope and SVS files, but the trick is importing them into Matlab.

In this case, we’ll run through a quick 2 region mask to give the general approach and code structure needed to create binary masks. Here we can see 2 green annotations made using the pen tool:

annotated

ImageScope saves the annotation x,y coordinates relative to the base image of the pyramid (as discussed here). It stores them in an xml file of the same name as the original file except with an xml extension. So the steps are rather straight forward from there

  1. Parse xml file
  2. Extract x,y coordinates
  3. Make binary mask using x,y points

One caveat, which we discussed in the previous tutorial, is that the base image is likely too large to be loaded into memory. This means that making a binary image of the same size as the binary image is also unlikely to work. To compensate we’ll scale down the binary mask to a different layer in the image pyramid, one which can/does fit easily into memory.

Looking at the XML

Lets see what the xml looks like:

xml

We can see some interesting pieces of information, such as the microns per pixel, but we’re more focused on the region annotations. Here we can see two regions, with ID=”1” and ID=”2” (lines 14 and 392 respectively). The attributes associated with the region contain things like length, area, etc. More importantly, we can see region 2 is made up of a set of vertices consisting of x,y coordinates. Our goal is to be able to extract these in an orderly fashion.

Extracting the XML

Since the xml is rather well structured, we can use the underlying xerces dom class to extract the data. The code looks something like this:

  1. svs_file='TCGA-A1-A0SD-01Z-00-DX1.DB17BFA9-D951-42A8-91D2-F4C2EBC6EB9F.svs';
  2. xml_file='TCGA-A1-A0SD-01Z-00-DX1.DB17BFA9-D951-42A8-91D2-F4C2EBC6EB9F.xml';
  3.  
  4.  
  5. xDoc = xmlread(xml_file);
  6. Regions=xDoc.getElementsByTagName('Region'); % get a list of all the region tags
  7. for regioni = 0:Regions.getLength-1
  8.     Region=Regions.item(regioni);  % for each region tag
  9.  
  10.     %get a list of all the vertexes (which are in order)
  11.     verticies=Region.getElementsByTagName('Vertex');
  12.     xy{regioni+1}=zeros(verticies.getLength-1,2); %allocate space for them
  13.     for vertexi = 0:verticies.getLength-1 %iterate through all verticies
  14.         %get the x value of that vertex
  15.         x=str2double(verticies.item(vertexi).getAttribute('X'));
  16.        
  17.         %get the y value of that vertex
  18.         y=str2double(verticies.item(vertexi).getAttribute('Y'));
  19.         xy{regioni+1}(vertexi+1,:)=[x,y]; % finally save them into the array
  20.     end
  21.    
  22. end

 

Afterwards, we can plot the results using code like this:

  1. figure,hold all
  2. set(gca,'YDir','reverse'); %invert y axis
  3. for zz=1:length(xy)
  4.     plot(xy{zz}(:,1),xy{zz}(:,2),'LineWidth',12)
  5. end

The only issue is that we need to invert the y-axis because images have the upper left corner as (0,0) and matlab’s default plotting uses the bottom left corner as the origin. The result of the plotting look like this, which are exactly the same shape as our original annotations (!) :

plot_xy_thick

 

 

Making the Mask

Now the issue is that we have a plot but not a binary mask. Maybe we can just use them directly, lets try to allocate the necessary memory:

max_xy

 

Okay that’s not going to work, we have the same issue as we’ve discussed before. The base image is simply too large to fit into memory, and while the binary mask is 1/3 the size (single channel versus RGB), its still too large for my 16GB of RAM.

 

So what if we still want to look at it? Well we can apply the same trick as before and scale the points to match a smaller pyramidal layer.

  1. svsinfo=imfinfo(svs_file);
  2. s=1; %base level of maximum resolution
  3. s2=5; % down sampling of 1:32
  4. hratio=svsinfo(s2).Height/svsinfo(s).Height;  %determine ratio
  5. wratio=svsinfo(s2).Width/svsinfo(s).Width;
  6.  
  7. nrow=svsinfo(s2).Height;
  8. ncol=svsinfo(s2).Width;
  9. mask=zeros(nrow,ncol); %pre-allocate a mask
  10. for zz=1:length(xy) %for each region
  11.     smaller_x=xy{zz}(:,1)*wratio; %down sample the region using the ratio
  12.     smaller_y=xy{zz}(:,2)*hratio;
  13.    
  14.     %make a mask and add it to the current mask
  15.     %this addition makes it obvious when more than 1 layer overlap each
  16.     %other, can be changed to simply an OR depending on application.    
  17.     mask=mask+poly2mask(smaller_x,smaller_y,nrow,ncol);
  18. end

From here we can simply display the mask with imshow(mask):

mask

And see that the interior is indeed the same as what was marked above in the ImageScope version.

All done!

Full source available here

7 thoughts on “Working with Aperio SVS files in Matlab – Converting Annotations to Binary Masks”

    1. I’ve never had access to such a machine so unfortunately I don’t know :-\ If you find out, I’d appreciate any insights you may discover!

  1. Have you used only the level s3=3 for svs spliting in all of your researches?

    Is it a way to make precise xml mask for base level s=1
    maybe by:

    1)Making a window sized 20*20 pixels.
    2)Iterate it through mask at resolution s2=5
    3)>mask2poly_for_window=transpose(bwboundaries(mask(window(i))))
    4)Rescale poly for a piece of svs, backwards to s=1
    multipling components of mask2poly_for_window to 1/hratio(..) and 1/wratio(…) correspondingly.

    We would have enough memory, because we save not a huge mask of all svs (9*10^4×9*10^4) but a small piece of it.
    The price is a big amount of iterations.

    Is there any need in such procedure or basically we should need more maxpools further if we use the same operating memory?

    1. we use all levels depending on the necessary magnifications.

      That said, for very large images (WSI), we’ve found that its far easier to extract smaller regions which are manageable in size, allowing for easier parallelization and viewing. Reading and writing individual tiles tends to be extremely slow, and error prone as it requires a number of wrappers to address the various formats/styles (openslide, etc). That is to say, even if I had a 100,000 x 100,000 mask for an associated 100k x 100k image, I’m not entirely sure what I would do with it :)

  2. Line 10 in the last code block, no need to use the loop. The poly2mask function does things in one go.

    Great work any way. Very helpful.

    1. i think you’re missing the case where there is more than 1 annotated region. poly2mask will *connect* all of the vertices in the list, thus merging together two regions. the for loop is to generate them separately and combined

Leave a Reply

Your email address will not be published. Required fields are marked *