<jupyter_start><jupyter_text>QC script to calculate sharpness metric for images in a plate The following script will access the IDR images ina facility manager's context.The QC script does the following:1. Extracts Images from IDR (read-only)2. Calculates sharpness based on the algorithm definedin the following publication  * [Image Sharpness Measure for Blurred Images in Frequency Domain](https://doi.org/10.1016/j.proeng.2013.09.086)3. Creates a numpy array of the sharpness scores for every well,  * separately for every field and every channel  * the numpy array is then reshaped to the same dimensions of the plate   e.g. 96 well plates will create an array with 8 rows and 12 columns.4. Plots a heatmap for every field and every channel,   and arranges all plots within a subplot.  Workflow summary ![Overview](./includes/PlateToHeatmap.jpg) ***Import Packages***<jupyter_code>import string
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from skimage import feature
from scipy import ndimage
from scipy.ndimage import convolve
from scipy import misc
from IPython.display import HTML
from idr import connection<jupyter_output><empty_output><jupyter_text>Read data from IDR (Public Database : Read-only!)  Create connection and plate identifier<jupyter_code>conn = connection('idr.openmicroscopy.org')
plateId = 408<jupyter_output>Connected to IDR ...<jupyter_text>Fetch plate object and print details<jupyter_code>plate = conn.getObject("Plate", plateId)
print("\nNumber of fields:", plate.getNumberOfFields())
print("\nGrid size:", plate.getGridSize())
print("\nWells in Plate:", plate.getName())

plate_rows = plate.getRows()
plate_columns = plate.getColumns()
plate_format = plate_rows * plate_columns

print("\nPlate Format:", plate_format)<jupyter_output>Number of fields: (0L, 0L)

Grid size: {'rows': 16L, 'columns': 24L}

Wells in Plate: 0001-03--2005-08-01

Plate Format: 384<jupyter_text>***Algorithm List***<jupyter_code>class AlgorithmList:
    def fourierBasedSharpnessMetric(self):
        fftimage = np.fft.fft2(plane)
        fftshift = np.fft.fftshift(fftimage)
        fftshift = np.absolute(fftshift)
        M = np.amax(fftshift)
        Th = (fftshift > M // float(1000)).sum()
        if 'image' in locals():
            size = float(image.getSizeX() * image.getSizeY())
            sharpness = Th // size
            return sharpness * 10000
        else:
            return Th

    def gradientBasedSharpnessMetric(self):
        gy, gx = np.gradient(plane)
        gnorm = np.sqrt(gx**2 + gy**2)
        sharpness = np.average(gnorm)
        return sharpness

    def edgeBasedSharpnessMetric(self):
        edges1 = feature.canny(plane, sigma=3)
        kernel = np.ones((3, 3))
        kernel[1, 1] = 0
        sharpness = convolve(edges1, kernel, mode="constant")
        sharpness = sharpness[edges1 != 0].sum()
        return sharpness


print("loaded:", dir(AlgorithmList))<jupyter_output>loaded: ['__doc__', '__module__', 'edgeBasedSharpnessMetric', 'fourierBasedSharpnessMetric', 'gradientBasedSharpnessMetric']<jupyter_text>Test your algorithm on example data<jupyter_code>resultArray = np.zeros((5, 2), dtype=float)
plt.figure(figsize=(20, 15))
cntr = 1
for sigValue in range(0, 20, 4):
    face = misc.face(gray=True)
    plane = ndimage.gaussian_filter(face, sigma=sigValue)
    plt.subplot(1, 5, cntr)
    plt.imshow(plane, cmap=plt.cm.gray)
    plt.axis('off')
    sharpness = AlgorithmList().fourierBasedSharpnessMetric()
    resultArray[cntr - 1, 1] = sharpness
    resultArray[cntr - 1, 0] = sigValue
    cntr = cntr + 1

plt.show()
plt.figure(figsize=(10, 8))
plt.plot(resultArray[:, 0], resultArray[:, 1], 'ro')
plt.xlabel('Levels of gaussian blur')
plt.ylabel('sharpness score')
plt.show()

plt.gcf().clear()
print(resultArray)<jupyter_output><empty_output><jupyter_text>Test your algorithm on plate data<jupyter_code>imageId = 171499
image = conn.getObject("Image", imageId)
print(image.getName(), image.getDescription())

pixels = image.getPrimaryPixels()
image_plane = pixels.getPlane(0, 0, 0)

resultArray = np.zeros((5, 2), dtype=float)
plt.figure(figsize=(20, 15))
cntr = 1
for sigValue in range(0, 20, 4):
    face = misc.face(gray=True)
    plane = ndimage.gaussian_filter(image_plane, sigma=sigValue)
    plt.subplot(1, 5, cntr)
    plt.imshow(plane, cmap=plt.cm.gray)
    plt.axis('off')
    sharpness = AlgorithmList().fourierBasedSharpnessMetric()
    resultArray[cntr - 1, 1] = sharpness
    resultArray[cntr - 1, 0] = sigValue
    cntr = cntr + 1
plt.show()
plt.figure(figsize=(10, 8))
plt.plot(resultArray[:, 0], resultArray[:, 1], 'ro')
plt.xlabel('Levels of gaussian blur')
plt.ylabel('sharpness score')
plt.show()
plt.gcf().clear()
print(resultArray)<jupyter_output>0001-03--2005-08-01 [Well 125, Field 1 (Spot 125)]<jupyter_text>Iterative calculations for the whole plate`stride_r = 4` and `stride_c = 4` so that this notebookcan be run quickly as an example.To check all 384 wells change them both to `1`.<jupyter_code>stride_r = 4
stride_c = 4
chnames = None
cntr = 0
fields = 0
size_z = fields
print("Iterating through wells...")
rc = dict(((well.row, well.column), well) for well in plate.listChildren())
for wrc in sorted(rc.keys()):
    if wrc[0] % stride_r > 0 or wrc[1] % stride_c > 0:
        continue
    well = rc[wrc]
    print('Row:%d Column:%d' % (well.row, well.column))
    index = well.countWellSample()
    image = well.getImage(fields)
    if chnames is None:
        chnames = [ch.getLabel() for ch in image.getChannels(True)]
    pixels = image.getPrimaryPixels()
    size_c = image.getSizeC()
    if cntr == 0:
        result_array = np.full((plate_format, size_c), np.nan)
    for ch in range(0, size_c):
        plane = pixels.getPlane(0, ch, 0)
        sharpness = AlgorithmList().fourierBasedSharpnessMetric()
        wellid = well.row * plate_columns + well.column
        result_array[wellid, ch] = sharpness
        tempvalue = result_array[wellid, ch]
        fieldid = (fields + ch * size_c)
        cntr = cntr + 1<jupyter_output>Iterating through wells...
Row:0 Column:0
Row:0 Column:4
Row:0 Column:8
Row:0 Column:12
Row:0 Column:16
Row:0 Column:20
Row:4 Column:0
Row:4 Column:4
Row:4 Column:8
Row:4 Column:12
Row:4 Column:16
Row:4 Column:20
Row:8 Column:0
Row:8 Column:4
Row:8 Column:8
Row:8 Column:12
Row:8 Column:16
Row:8 Column:20
Row:12 Column:0
Row:12 Column:4
Row:12 Column:8
Row:12 Column:12
Row:12 Column:16
Row:12 Column:20<jupyter_text>Reshape numpy array and plot heatmaps<jupyter_code>alphabets = list(string.ascii_uppercase)
plate_name = plate.getName()
colval = 0
planes = []
cntr = 0
size_c = 3
fig = plt.figure(figsize=(30, 15))
for rowval in range(0, size_c):
    data = result_array[:, rowval].reshape(plate_rows, plate_columns)
    ax = plt.subplot(size_c, 1, cntr + 1)
    plt.pcolor(data)
    plt.colorbar()
    ax.title.set_text(chnames[rowval])
    plt.xticks(np.arange(0.5, plate_columns, 1.0))
    plt.yticks(np.arange(0.5, plate_rows, 1.0))
    xlabels = range(1, plate_columns+1)
    ax.set_xticklabels(xlabels)
    ylabels = range(1, plate_rows+1)
    ax.set_yticklabels([alphabets[i - 1] for i in ylabels])
    plt.gca().invert_yaxis()
    plt.clim(0, 40000)

    data = np.repeat(data, 20, axis=1)
    data = np.repeat(data, 20, axis=0)
    planes.append(np.uint16(data))
    cntr = cntr + 1

plt.show()
fig.savefig(plate_name + 'SharpnessHeatMaps.png')<jupyter_output><empty_output><jupyter_text>***Thumbnails of top2 and bottom 2 percentile images***<jupyter_code>mapAnnotationNameSpace = "openmicroscopy.org/mapr/gene"
bulkAnnotationNameSpace = "openmicroscopy.org/omero/bulk_annotations"


def id_to_image_html(id):
    return '<img src="http://idr.openmicroscopy.org/webclient/render_thumbnail/%d/"/>' % id


def getGeneInformation(image):
    id = image.getId()
    image1 = conn.getObject('Image', id)
    cc = image1.getAnnotation(mapAnnotationNameSpace)
    rows = cc.getValue()
    html = []
    for r in rows:
        if r[1].startswith("http"):
            tempvar = "<a href='" + r[1] + "'>" + r[1] + "</a>"
        else:
            tempvar = r[1]
        html.append("<tr><td>" + tempvar + "</td></tr>")
    return ("<table>" + "".join(html) + "</table>")


def getQualityControl(image):
    id = image.getId()
    image1 = conn.getObject('Image', id)
    cc = image1.getAnnotation(bulkAnnotationNameSpace)
    rows = cc.getValue()
    html = []
    for r in rows:
        if r[0].startswith('Control') or r[0].startswith('Quality'):
            html.append("<tr><td>" + r[1] + "</td></tr>")
    return ("<table>" + "".join(html) + "</table>")


fields = 0
ch = 2

result_array_ch = result_array[:, ch]
threshold = np.percentile(result_array_ch[~np.isnan(result_array_ch)], 2)
imageList = []
for wrc in sorted(rc.keys()):
    if wrc[0] % stride_r > 0 or wrc[1] % stride_c > 0:
        continue
    well = rc[wrc]
    row = well.row
    column = well.column
    sharpness = result_array[((row)*plate_columns) + column, ch]
    if (sharpness <= threshold):
        image = well.getImage(fields)
        imageList.append(image)

images = [(x.id, x.id, x.getName(), x, x) for x in (imageList)]
pd.set_option('display.max_colwidth', -1)
df = pd.DataFrame(images, columns=['Id', 'Image', 'Name',
                                   'GeneInformation',
                                   'QualityControl'])
HTML(df.to_html(escape=False,
                formatters=dict(Image=id_to_image_html,
                                GeneInformation=getGeneInformation,
                                QualityControl=getQualityControl)))
mapAnnotationNameSpace = "openmicroscopy.org/mapr/gene"
bulkAnnotationNameSpace = "openmicroscopy.org/omero/bulk_annotations"


def id_to_image_html(id):
    return '<img src="http://idr.openmicroscopy.org/webclient/render_thumbnail/%d/"/>' % id


def getGeneInformation(image):
    id = image.getId()
    image1 = conn.getObject('Image', id)
    cc = image1.getAnnotation(mapAnnotationNameSpace)
    rows = cc.getValue()
    html = []
    for r in rows:
        if r[1].startswith("http"):
            tempvar = "<a href='" + r[1] + "'>" + r[1] + "</a>"
        else:
            tempvar = r[1]
        html.append("<tr><td>" + tempvar + "</td></tr>")
    return ("<table>" + "".join(html) + "</table>")


def getQualityControl(image):
    id = image.getId()
    image1 = conn.getObject('Image', id)
    cc = image1.getAnnotation(bulkAnnotationNameSpace)
    rows = cc.getValue()
    html = []
    for r in rows:
        if r[0].startswith('Control') or r[0].startswith('Quality'):
            html.append("<tr><td>" + r[1] + "</td></tr>")
    return ("<table>" + "".join(html) + "</table>")


fields = 0
ch = 2
result_array_ch = result_array[:, ch]
threshold = np.percentile(result_array_ch[~np.isnan(result_array_ch)], 98)
imageList = []
for well in plate.listChildren():
    row = well.row
    column = well.column
    sharpness = result_array[((row)*plate_columns) + column, ch]
    if (sharpness >= threshold):
        image = well.getImage(fields)
        imageList.append(image)

images = [(x.id, x.id, x.getName(), x, x) for x in (imageList)]
pd.set_option('display.max_colwidth', -1)
df = pd.DataFrame(images, columns=['Id', 'Image',
                                   'Name', 'GeneInformation',
                                   'QualityControl'])
HTML(df.to_html(escape=False,
                formatters=dict(Image=id_to_image_html,
                                GeneInformation=getGeneInformation,
                                QualityControl=getQualityControl)))<jupyter_output><empty_output><jupyter_text>Close Connection to IDR<jupyter_code>conn.close()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1387298059663219, 'fraction_numerical': 0.026327036498846056, 'mean_word_length': 3.582843713278496, 'pattern_counts': {'":': 0, '<': 67, '<?xml version=': 0, '>': 72, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '414609', 'license': 'BSD-2-Clause', 'n_tokens_mistral': 4265, 'n_tokens_neox': 3980, 'n_words': 935, 'path': 'CalculateSharpness.ipynb', 'repo_name': 'IDR/idr-notebooks', 'size': '11699'}