volume-visualization/imageprocessing/sliceMapCreator.py

import io
from skimage import img_as_ubyte
from skimage.transform import resize

from PIL import Image

from django.conf import settings

from multiprocessing import Process, cpu_count, Queue

import numpy as np
import math
import copy

import logging

import pdb

logger = logging.getLogger(__name__)


class Sprite(object):
    imagesAsArrays = None
    imagesToProcess = None
    imagesPerRow = 0
    currentRow = 0
    imagesInCurRow = 0
    totalRows = 0
    numberOfImages = 0
    numberOfImagesInSprite = 0
    destFilename = None
    buffers = []
    contenttype = None,
    numberOfSprites = 0
    slicesPerSprite = 0

    thresholdIndex_otsu = 255
    thresholdIndex_yen  =  255
    thresholdIndex_isodata  =  255

def __initializeSprite(imageFrom, imageTo, sliceShape, dtype):

    imagesPerRow, totalRows, imagesToProcess, numberOfSprites, slicesPerSprite, zDownScalingFactor = __calculateRowsAndColsOfSpriteMap(imageFrom, imageTo)

    logger.debug('sliceMapCreator opened source image, writing all data into one array now')

    outputImageArrays = []
    for counter in range(0, numberOfSprites):
        outputImageArray = np.empty((sliceShape[0] * totalRows, sliceShape[1] * imagesPerRow), dtype=dtype)
        outputImageArray.setflags(write=True)
        outputImageArrays.append(outputImageArray)

    sprite = Sprite()
    sprite.imagesToProcess = imagesToProcess
    sprite.imagesPerRow = imagesPerRow
    sprite.totalRows = totalRows
    sprite.numberOfSprites = numberOfSprites
    sprite.slicesPerSprite = slicesPerSprite
    sprite.imagesAsArrays = outputImageArrays

    return sprite


def __calculateRowsAndColsOfSpriteMap(imageFrom, imageTo):
    imagesToProcess = imageTo - imageFrom + 1
    slicesPerSprite = imagesToProcess

    # upperbound are 25 ^ 2 slices per map
    maxSlicesPerMap = settings.MAX_NUMBER_SLICES_PER_MAP

    zDownScalingActivated = settings.Z_DOWNSCALING_ACTIVATED
    zDownScalingFactor = 1

    if zDownScalingActivated is False:
        # calculate the number of sprites for one frame
        numberOfSprites = int(math.ceil(imagesToProcess / maxSlicesPerMap))
    else:
        numberOfSprites = 1
        zDownScalingFactor = int(math.ceil(imagesToProcess / maxSlicesPerMap))

    imagesToProcess = int(math.floor(imagesToProcess / float(zDownScalingFactor)))

    # stretch the whole number of images on to the sprites
    slicesPerSprite = imagesToProcess / numberOfSprites

    imagesPerRow = int(math.ceil(math.sqrt(slicesPerSprite)))
    # now i has a squared number of images
    slicesPerSprite = math.pow(imagesPerRow, 2)
    maxImageCount = slicesPerSprite * numberOfSprites

    imagesLeft = maxImageCount - imagesToProcess
    logger.debug('sliceMapCreator: imagesToProcess: %d, emptyImages in the last rows: %d ' % (imagesToProcess, imagesLeft))

    totalRows = int(imagesPerRow)
    logger.debug('sliceMapCreator imagesPerRow: %d, totalRows: %d' % (imagesPerRow, totalRows))

    return (imagesPerRow, totalRows, imagesToProcess, numberOfSprites, slicesPerSprite, zDownScalingFactor)


def __copyImageToSprite(sliceImage, destArray, sprite):
    numY, numX = sliceImage.shape

    pixelRow = sprite.imagesInCurRow * numX
    pixelCol = sprite.currentRow * numY

    destArray[pixelCol:pixelCol+numY, pixelRow:pixelRow+numX] = sliceImage

    return destArray


def processFrames(frames, frameFrom, frameTo, imageFrom, imageTo, imgFormat, destSize):

    if len(frames.shape) != 4:
        raise ValueError('due to memory restrictions and new implementation, currently only frames of 4 dimensions supported')

    sprites = []
    imagesPerRow, totalRows, imagesToProcess, numberOfSprites, slicesPerSprite, zDownScalingFactor = __calculateRowsAndColsOfSpriteMap(imageFrom, imageTo)

    upperBoundPixels = 16384

    if frames.shape[3] * imagesPerRow > upperBoundPixels:
        pixelsPerImage = int(math.ceil(upperBoundPixels / imagesPerRow))
        newImageShape = (pixelsPerImage, pixelsPerImage)
        newFrames = np.empty((frames.shape[0], frames.shape[1], newImageShape[0], newImageShape[1]), dtype=frames.dtype)
        frameCounter = 0

        logger.debug('sliceMapCreator resizing slices to (%d, %d)' % newImageShape)
        for frameCounter in range(0, frames.shape[0]):
            imageCounter = 0
            cpus = cpu_count()
            imagesPerCpu = int(math.ceil(frames.shape[1] / float(cpus)))
            processes = []
            queues = [Queue() for i in range(0, cpus)]

            for cpu in range(0, cpus):
                fromImage = imagesPerCpu * cpu
                toImage = (imagesPerCpu * (cpu + 1)) - 1
                if fromImage >= frames.shape[1]:
                    break

                if toImage >= frames.shape[1]:
                    toImage = frames.shape[1] - 1

                process = Process(target=framePreprocessor, args=(frames[frameCounter], queues[cpu], fromImage, toImage, newImageShape,))
                processes.append(process)

            logger.debug('sliceMapCreator before start loop')
            for process in processes:
                process.start()

            resultFrames = []
            logger.debug('sliceMapCreator empty queues')
            for queue in queues:
                resultFrames.append(queue.get())

            logger.debug('sliceMapCreator merging results')

            for resultFrame in resultFrames:

                imageCounter = resultFrame.fromImage
                frame = resultFrame.frame
                for image in frame:
                    newFrames[frameCounter][imageCounter] = image
                    imageCounter += 1

            for process in processes:
                process.join()

            '''
            old sequential implementation

            for imageCounter in range(0, frames.shape[1]):

                image = resize(frames[frameCounter][imageCounter], (pixelsPerImage, pixelsPerImage))
                newFrames[frameCounter][imageCounter] = img_as_ubyte(image)

                frames[frameCounter][imageCounter][:][:] = np.empty((frames.shape[2], frames.shape[3]), dtype=frames.dtype)[:][:]
            '''

        del frames
        frames = newFrames
    else:
        logger.debug('sliceMapCreator no need to resize images to lower shape')

    if frameFrom is not None and len(frames.shape) > 3:
        frameNumber = frameFrom
        for frame in frames:
            sprite = processFrame(frame, frameNumber, imageFrom, imageTo, zDownScalingFactor)
            sprites.append(sprite)
            del frame
            frameNumber += 1
    else:
        sprite = processFrame(frames, None, imageFrom, imageTo)
        sprites.append(sprite)

    del frames

    return __resizeSprites(sprites, destSize, imgFormat)


def processFrame(frame, frameNumber, imageFrom, imageTo, zDownScalingFactor):
    logger.debug('sliceMapCreator frameNumber %d' % frameNumber)

    if len(frame) is 0:
        raise ValueError('could not find any frame')

    firstSlice = frame[0]

    sprite = __initializeSprite(imageFrom, imageTo, firstSlice.shape, frame.dtype)

    imageCounter = 0
    frameLength = len(frame)

    while imageCounter < frameLength:
        image = frame[imageCounter]

        arrayIndex = int(math.floor(float(sprite.numberOfImages) / sprite.slicesPerSprite))

        __copyImageToSprite(image, sprite.imagesAsArrays[arrayIndex], sprite)

        sprite.imagesInCurRow += 1
        sprite.numberOfImages += 1
        sprite.numberOfImagesInSprite += 1

        imageCounter += (1 * zDownScalingFactor)
        if sprite.imagesInCurRow == sprite.imagesPerRow:
            if sprite.numberOfImages >= sprite.imagesToProcess:
                logger.debug('sliceMapCreator copy end reached, breaking now')
                break

            sprite.currentRow += 1
            sprite.imagesInCurRow = 0

            if sprite.numberOfImagesInSprite == sprite.slicesPerSprite:
                logger.debug('sliceMapCreator resetting currentrow to 0')
                sprite.currentRow = 0
                sprite.numberOfImagesInSprite = 0

    sprite = calculateThresholdIndexes(sprite)

    return sprite

def calculateThresholdIndexes(sprite):
    from skimage.filters import threshold_otsu, threshold_yen, threshold_isodata

    # pdb.set_trace()
    
    image = sprite.imagesAsArrays[0]

    sliceMapWidth = len(image[0])
    sliceMapHeight = len(image)

    numberSlicesInSliceMap = sprite.numberOfImages if sprite.numberOfImages >= 1 else 1

    RowInSliceMap = sprite.imagesPerRow
    ColInSliceMap = sprite.imagesPerRow

    sliceWidth = sliceMapWidth / ColInSliceMap;
    sliceHeight = sliceMapHeight / RowInSliceMap;

    numberOfCentralSlice = (numberSlicesInSliceMap - 1) / 2

    x = numberOfCentralSlice / RowInSliceMap
    xp = x * sliceWidth

    y = numberOfCentralSlice - (x * RowInSliceMap)
    yp = y * sliceHeight

    centralSliceX0 = xp
    centralSliceX1 = xp + sliceWidth

    centralSliceY0 = yp
    centralSliceY1 = yp + sliceHeight

    centerSliceImage = image[centralSliceX0 : centralSliceX1, centralSliceY0 : centralSliceY1]

    sprite.thresholdIndex_otsu = threshold_otsu(centerSliceImage)
    sprite.thresholdIndex_yen  =  threshold_yen(centerSliceImage)
    sprite.thresholdIndex_isodata  =  threshold_isodata(centerSliceImage)

    print("thresholdIndex_otsu = " + str(sprite.thresholdIndex_otsu))
    print("thresholdIndex_yen = " + str(sprite.thresholdIndex_yen))
    print("thresholdIndex_isodata = " + str(sprite.thresholdIndex_isodata))

    return sprite

def generateInfo(sprites, frameFrom, frameTo, imageFrom, imageTo):
    logger.debug('sliceMapCreator generating and saving .js info file')

    firstSprite = list(sprites.values())[0][0]
    infoObject = {}
    infoObject['slicesOverX'] = firstSprite.imagesPerRow
    infoObject['slicesOverY'] = firstSprite.totalRows
    infoObject['numberOfSlices'] = firstSprite.numberOfImages
    infoObject['imgMin'] = imageFrom
    infoObject['imgMax'] = imageTo
    infoObject['numberOfSprites'] = firstSprite.numberOfSprites
    infoObject['slicesPerSprite'] = firstSprite.slicesPerSprite
    multipleFrames = len(sprites) > 1
    infoObject['multipleFrames'] = multipleFrames
    if multipleFrames:
        infoObject['frameFrom'] = frameFrom
        infoObject['frameTo'] = frameTo

    infoObject['thresholdIndex_otsu'] = firstSprite.thresholdIndex_otsu
    infoObject['thresholdIndex_yen'] = firstSprite.thresholdIndex_yen
    infoObject['thresholdIndex_isodata'] = firstSprite.thresholdIndex_isodata

    return infoObject


def __resizeSprites(sprites, destSize, imgFormat):
    returnSprites = {}
    exponentMin = 9
    exponent = exponentMin

    firstSprite = sprites[0]
    # hardcoded value, 2^12
    # 2^13 webp files force error during saving
    widthUpperBound = 4096

    generatorSprite = copy.deepcopy(firstSprite)
    del generatorSprite.imagesAsArrays

    logger.debug('sliceMapCreator resizeSprites, shape of input sprites is (%d, %d)' % sprites[0].imagesAsArrays[0].shape)

    if destSize is not None:
        widthUpperBound = destSize

    exponentFound = False

    # calculate the power of 2 size of the sprite map
    while not exponentFound:
        pot = math.pow(2, exponent)
        if pot > widthUpperBound:
            exponentFound = True
            exponent = exponent
        else:
            exponent += 1

    for exponentToStore in range(exponent - 1, exponentMin - 1, -1):
        destSize = math.pow(2, exponentToStore)
        logger.debug('sliceMapCreator resize sprites to %d' % destSize)

        for sprite in sprites:
            buffers = []
            for arrayCounter in range(0, sprite.numberOfSprites):
                imageArray = sprite.imagesAsArrays[arrayCounter]
                # for lower memory consumption, we resize the image continuesly --> written right?
                imageArray = resize(imageArray, (destSize, destSize))
                imageArray = img_as_ubyte(imageArray)

                buff = io.BytesIO()

                quality = 95
                imgFormat = imgFormat.lower()
                pilImage = Image.fromarray(imageArray)

                contenttype = 'image/'

                if imgFormat == 'webp':
                    pilImage = pilImage.convert('RGBA')
                    try:
                        pilImage.save(buff, format='WEBP', quality=quality)
                    except OSError:
                        # due to a bug in the library of webp, the webp encoder fails for high file sizes
                        # --> if this is a big 'destSize', we have to try it again with lower quality levels
                        logger.debug('sliceMapCreator failed to store webp with size of %dpx and quality value of %d' % (destSize, quality))
                        logger.debug('sliceMapCreator will not use it, skipping further processing of %dpx sprite' % destSize)

                        continue
                    contenttype += 'webp'

                elif imgFormat == 'png':
                    pilImage.save(buff, format='PNG')
                    contenttype += 'png'

                else:
                    pilImage.save(buff, format='JPEG', quality=quality, progressive=True)
                    contenttype += 'jpeg'

                logger.debug('sliceMapCreator image with res %d has size of: %d' % (destSize, buff.tell()))

                buff.seek(0, 0)
                buffers.append(buff)

            # end for each array

            spriteToAdd = copy.deepcopy(generatorSprite)
            spriteToAdd.buffers = buffers
            spriteToAdd.contenttype = contenttype

            # sprite.image
            if '%d' % destSize not in returnSprites:
                returnSprites['%d' % destSize] = []

            returnSprites['%d' % destSize].append(spriteToAdd)

    for sprite in sprites:
        del sprite.imagesAsArrays

    return returnSprites


def calculateXYDimensions(frame):
    if len(frame.shape) != 3:
        raise ValueError('frame must have the dimension 3')

    new_height = frame.shape[1]
    new_width = frame.shape[2]

    max_size = max(new_height, new_width)

    pixel_offset = [0, 0]

    if new_width < max_size:
        pixel_offset[0] = math.ceil((max_size - new_width) / 2)
        new_width = max_size

    if new_height < max_size:
        pixel_offset[1] = math.ceil((max_size - new_height) / 2)
        new_height = max_size

    logger.debug('sliceMapCreator pixeloffsets %d, %d, newWidth %d' % (pixel_offset[0], pixel_offset[1], new_width))

    return pixel_offset, (new_width, new_height)


def crop_frames(frames, pixel_offset, new_size):
    if len(frames.shape) != 4:
        raise ValueError('frame must have the dimension 3')

    offset_x = pixel_offset[0]
    offset_y = pixel_offset[1]

    t, z, y, x = frames.shape

    squared_frames = np.empty((t, z, new_size[1], new_size[0]), dtype=frames.dtype)

    logger.debug('squared_frames.shape: %d, %d, %d, %d' % squared_frames.shape)
    logger.debug('        frames.shape: %d, %d, %d, %d' % frames.shape)
    logger.debug('offsety: %d, offsetx: %d ' % (offset_y, offset_x))
    '''for frameNumber in range(0, frames.shape[0]):
        squared_frames[frameNumber] = sliceMapCreator.cropFrame(frames[frameNumber], pixelOffset, newSize)'''

    squared_frames[:, :, offset_y: y + offset_y, offset_x: x + offset_x] = frames

    '''newFrame = np.zeros((z, y + 2 * offset_y, x + 2 * offset_x), dtype=frame.dtype)
    newFrame[:, offset_y : offset_y + y, offset_x : offset_x + x]'''

    return squared_frames


class ResultFrame():
    def __init__(self, frame, fromImage):
        self.frame = frame
        self.fromImage = fromImage


def framePreprocessor(inFrame, queue, fromImage, toImage, destShape):
    numberOfImages = toImage - fromImage + 1
    outFrame = np.empty((numberOfImages, destShape[0], destShape[1]), dtype=inFrame.dtype)
    outputCounter = 0
    for imageNumber in range(int(fromImage), int(toImage) + 1):
        image = resize(inFrame[imageNumber], destShape)
        outFrame[outputCounter] = img_as_ubyte(image)
        outputCounter += 1

    queue.put(ResultFrame(outFrame, fromImage))