In this post
@kitcat published virtual slide images that contain checkerboard pattern structure on the high-resolution planes.
Instead of speculating about the meaning of this patterns I want to show here that they can be reduced or even eliminated by Gaussian smoothing while the image content and the de facto resolution is not reduced significantly.
I have used @petebankhead 's great QuPath to save smoothed regions of the images in full resolution.
This approach seems to be appropriate because the image pixel size (0.25 µm, highest resolution) is less than the Abbe limit (~0.35 µm for 40x) and therefore the sampling is close to the Nyquist condition. The approach described here can be interpreted as kind of ‘reconstruction by smoothing’.
The resulting images show that the patterns do not really interfere with the visual use of the images, as they can be easily removed without essential loss of information by smooting with clij smoothing parameter = 0.7.
RGB color display
And the resulting images show that removing these patterns dramatically improves the quality of the images for digital processing by smooting with clij smoothing parameter = 1.0.
Red color display
DAB absorbance display
The following QuPath script
(saveSmoothedImageRegions.txt (11.6 KB) download file and change file extension to .groovy)
utilizes @haesleinhuepf ‘s clij/clupath GPU support
to save smoothed versions of image regions defined by annotations.
The script is designed for the processing of transmitted light color RGB images but it can be modified to work with single channel fluorescence images.
More information can be found in the script file.
For all RGB images loaded in the project, the script saves a pyramidal OME tiff image for each region defined in the image.
USAGE:
- Create an empty QuPath project
- Add relevant images to the project
- Select (multiple) annotations in each of the images
- Change to name of the annotations according to your schema (i.e. ‘sectionA1’, ‘sectionA2’, ‘sectionB1’, ‘sectionB2’ …)
- Load the script in the QuPath script editor
- Configure the parameters ‘outpath’, ‘outputDownsample’, ‘gaussSigma’ in the section PARAMETERS in the script
- Run the script in the QuPath script editor
- Wait until finished
Question @haesleinhuepf
Currently the data is first pulled into in ImagePlus and from there to the pixels array.
...
imp = clijx.pull(result)
pixels = (float[]) imp.getProcessor().getPixels()
...
Is there a more direct and faster way to pull the result directly to the pixels array?
Here is the relevant part of the script:
// ****** PARAMETERS *******
def outpath = "C:\\tmp\\QP_smoothed_out\\"
// OME.Tiff parameters
def tilesize = 512
def outputDownsample = 1
def pyramidscaling = 2
def compression = OMEPyramidWriter.CompressionType.J2K_LOSSY //J2K //UNCOMPRESSED //LZW
// CLIJ gauss smoothing parameter
double gaussSigma = 1.0
// ****** SCRIPT *******
Project project = QP.getProject()
for (entry in project.getImageList()) {
def name = entry.getImageName()
// Strip of file extension
name = name.substring(0, name.length() - 4)
print name
def ImageData currentImage = entry.readImageData()
def ImageServer currentServer = currentImage.getServer()
Collection<PathObject> annos = currentImage.getHierarchy().getAnnotationObjects()
for (anno in annos) {
def annoname = anno.getName()
annoRoi = anno.getROI()
int x = annoRoi.getBoundsX()
int y = annoRoi.getBoundsY()
int w = annoRoi.getBoundsWidth()
int h = annoRoi.getBoundsHeight()
ImageRegion selection = new ImageRegion(x, y, w, h, 1, 1);
def serverRGBsmooth = new TypeConvertServerRGBsmooth(currentServer, gaussSigma)
def transformedServer = new TransformedServerBuilder(serverRGBsmooth)
def myServer = transformedServer
.crop(selection)
//.rotate(RotatedImageServer.Rotation.ROTATE_90)
.build()
def filename = name + "_" + annoname
filename = filename + "_" + gaussSigma + "_" + outputDownsample+ "_" + tilesize + ".ome.tif"
//filename = filename + ".ome.tif"
println 'Writing OME-TIFF ' + filename
def pathOutput = outpath + filename
new OMEPyramidWriter.Builder(myServer)
.compression(compression)
.parallelize()
.tileSize(tilesize)
.scaledDownsampling(outputDownsample, pyramidscaling)
.build()
.writePyramid(pathOutput)
println('Done:' + filename)
}
}
println('Finished!')
// ****** CLASS DEFINITION : Smoothing RGB Server *******
class TypeConvertServerRGBsmooth extends TransformingImageServer<BufferedImage> {
ImagePlus imp
CLUPATH clijx
double gaussSigma
ColorModel cm
ImageServer<BufferedImage> currentserver
ArrayList<ColorTransformer.ColorTransformMethod> colortransformation = new ArrayList<ColorTransformer.ColorTransformMethod>()
TypeConvertServerRGBsmooth(ImageServer<BufferedImage> server, double gaussSigma) {
super(server)
currentserver = server
this.gaussSigma = gaussSigma
colortransformation.add(ColorTransformer.ColorTransformMethod.Red)
colortransformation.add(ColorTransformer.ColorTransformMethod.Green)
colortransformation.add(ColorTransformer.ColorTransformMethod.Blue)
cm = ColorModelFactory.getProbabilityColorModel32Bit(currentserver.getMetadata().getChannels())
// CLIJX
// Request 1st GPU
clijx = CLUPATH.getInstance()
// or request a specific GPU
// clijx = CLUPATH.getInstance("GeForce")
print(clijx.getGPUName() + '\n')
}
public ImageServerMetadata getOriginalMetadata() {
return currentserver.getMetadata()
}
@Override
protected ImageServerBuilder.ServerBuilder<BufferedImage> createServerBuilder() {
return currentserver.builder()
}
@Override
protected String createID() {
return UUID.randomUUID().toString();
}
@Override
public String getServerType() {
return "My RGB Type smoothing image server"
}
public BufferedImage readBufferedImage(RegionRequest request) throws IOException {
int idx, idxRoi
String path = request.getPath()
int ds = request.getDownsample()
int z = request.getZ()
int t = request.getT()
// Number of pixels of region extension
int extNP = Math.max(3, (int)(3 * gaussSigma))
int extNPminX, extNPmaxX, extNPminY, extNPmaxY
extNPminX = extNPmaxX = extNPminY = extNPmaxY = extNP
int xExt = request.x - extNPminX * ds
int yExt = request.y - extNPminY * ds
if (xExt < 0) {
extNPminX = (int) (1.0 * request.x / ds)
xExt = request.x - extNPminX * ds
}
if (yExt < 0) {
extNPminY = (int) (1.0 * request.y / ds)
yExt = request.y - extNPminY * ds
}
int wExt = request.width + (extNPminX + extNPmaxX) * ds
int hExt = request.height + (extNPminY + extNPmaxY) * ds
//print 'request: ' + ds + ' ' + request.x + ' ' + request.y + ' ' + request.width + ' ' + request.height + '\n'
//print 'requestExt: ' + ds + ' ' + xExt + ' ' + yExt + ' ' + wExt + ' ' + hExt + '\n'
// Create extended region to avoid edge effects of gaussian scmoothing
ImageRegion region = ImageRegion.createInstance(xExt, yExt, wExt, hExt, z, t);
RegionRequest requestExt = RegionRequest.createInstance(path, ds, region);
def img = getWrappedServer().readBufferedImage(requestExt)
def raster = img.getRaster()
int nBands = raster.getNumBands()
int w = img.getWidth()
int h = img.getHeight()
int wRoi = w - extNPminX - extNPmaxX
int hRoi = h - extNPminY - extNPmaxY
SampleModel model = new BandedSampleModel(DataBuffer.TYPE_BYTE, wRoi, hRoi, nBands)
byte[][] bytes = new byte[nBands][wRoi*hRoi]
DataBufferByte buffer = new DataBufferByte(bytes, wRoi*hRoi)
WritableRaster raster2 = Raster.createWritableRaster(model, buffer, null)
float[] pixels, pixelsRoi
int[] rgb = img.getRGB(0, 0, w, h, null, 0, w)
for (int b = 0; b < nBands; b++) {
pixels = ColorTransformer.getSimpleTransformedPixels(rgb, colortransformation.get(b), null);
// CLIJX
ClearCLBuffer input = clijx.pushArray(pixels, w, h, 1)
ClearCLBuffer result = clijx.create(input)
clijx.gaussianBlur(input, result, gaussSigma, gaussSigma)
imp = clijx.pull(result)
pixels = (float[]) imp.getProcessor().getPixels()
// Crop original region from extended region
pixelsRoi = new float[wRoi * hRoi]
for (int y = extNPminY; y < h - extNPmaxY; y++) {
idx = y * w + extNPminX
idxRoi = (y - extNPminY) * wRoi
System.arraycopy(pixels, idx, pixelsRoi, idxRoi, wRoi)
}
raster2.setSamples(0, 0, wRoi, hRoi, b, pixelsRoi)
input.close()
result.close()
}
return new BufferedImage(cm, Raster.createWritableRaster(model, buffer, null), false, null)
}
}