add auto-tiling and fix weighting after merging

[renerichter]

I just briefly checked into the mechanism of the toolbox and added some-functions which might do the intended trick for:

• automatic calculation of necessary padding (and shape)
• re-weighting after merging

. The code is implemented for 3D-images (hence for a 2D-image eg [1,y,x] could be provided) and needs to be slightly adjusted to work for 2D or ND. Things to be done:

• proper docstrings
• catching user-input-errors
• deleting unnecessary functions and clean-up

Further points are mentioned in the opened issue. :)

[the-lay]

Great job! I very sorry I am slow to respond, but I promise to thoroughly go through all your changes this week.

Thank you for your work!

[renerichter]

Thank you very much and take your time. I think there are still several errors in there as well as points to be discussed. I just added a function for automatic calculation of a minimal overlap (if no overlap is given), but it is not yet implemented in the init of Tiler class. Further, I fixed some small errors in the merger function.

[the-lay]

I was moving to a new apartment. Sorry you had to wait for so long, but I've finally had a chance to look it over! Instead of communicating in both #6 and #7, let's continue here, since we are discussing the PR after all :)

In addition to making sure it's ndimensional/docstring/edge cases/tests/cleanups/etc, some more discussion points:

1. I think instead of calling Tiler.calculate_padding during merging we should save original shape as a property of Tiler. That way the code can also be a bit cleaner since Tiler.calculate_padding and Tiler.fix_data_shape do not really need passed elements. Moreover, they can be converted into protected methods (I don't think end users want/need to be aware of these methods) or not even a method but just put it directly into Tiler.__init__.

2. For the reweighting I interpreted the calculation steps to be like if one would calculate a center-of-mass kind of thing. Hence,

   • multiplying each tile with the window
   • summing the tiles into the final image
   • dividing by the summed weights window-function

   could reproduce the wanted normalization of the output. On the other hand, if the window-functions are normed and symmetric than the final division by the weights would just undo the prior weighting of the tiles. Hence I believe I have a thinking error somehow here. What do you think?

I was thinking about this last two evenings and to be honest I think I only got more confused :D Do we need to divide by the weights? Yes, the output will not be the same as the input, but that's the whole point of applying window function, isn't it. Maybe just division by Tiler.data_visits is more correct... I will try out different ways in the next few days and will update you then.

3. Normalization for flat-top does not work because it has negative values, see plot. We can remove it or otherwise the weights[weights < atol] = 1 mechanism in Merger.norm_by_weights should change.

Again, thanks for your work, glad to see that it's not just me who wants something like that :) Would you mind to push the last version of what you have right now?

Also, if you are busy, I'd be willing to continue this myself. Let me know!

[renerichter]

No worries. I was quite busy last week as well and hence kept you waiting with my answer as well. I think I pushed my latest version of the box already, so there is right now no improvement yet. To your points:

I think instead of calling Tiler.calculate_padding during merging we should save original shape as a property of Tiler. That way the code can also be a bit cleaner since Tiler.calculate_padding and Tiler.fix_data_shape do not really need passed elements. Moreover, they can be converted into protected methods (I don't think end users want/need to be aware of these methods) or not even a method but just put it directly into Tiler.init.

Totally agreed. For all functions I newly implemented I wanted to make sure that it can be used as more general functions for different tasks as well and hence I gave them the option to pass elements. But you are right. From viewpoint of the toolbox these are not necessary and typically users would not want to use these functions in different contexts. So let's do it as you suggest.

I was thinking about this last two evenings and to be honest I think I only got more confused :D Do we need to divide by the weights? Yes, the output will not be the same as the input, but that's the whole point of applying window function, isn't it. Maybe just division by Tiler.data_visits is more correct... I will try out different ways in the next few days and will update you then.

That's actually the salt in the soup. What we want to achieve is the following:

1. tile image eg in two tiles M1 and M2 with overlap region dx
2. multiply weighting window W(x) on both tiles with the weights of tile 1 (w1) and tile 2 (w2) respectively
3. (optional: process images)
4. merge images by making sure that sum of weights in overlap-region is always 1, hence: W(x)=w1(x)+w2(x)=1

If you just normalize by the amount of visits step 4 is only true if each weight already satisfies wi(x)=1. But this seems to be only the case for boxcar. Hence I suggested: Normalization by the total-sum of the weights W(x) applied to a (ndimensional) overlap region should assure the normalization W(x)=1. Why would one like to have this normalization? Imagine that no processing is done to the image and hence splitting + merging should perfectly reconstruct the original image, no matter what the normalization function is (hence could be asymmetric, even though basically all of the functions we use are symmetric). Hence: the overlap of the various tiles should at max recreate the original image and hence M(x)=w1(x)*M1(x)+w2(x)*M2(x) for x in dx. Thus if M1 and M2 are just the same (hence: M1=M2=M ) and non-zero, then we could say: M(x)/M1(x)=1=w1(x)+w2(x). Hence: division by the weight map should be the right solution. what do you think?

I am back on normal track this week and would like to support you integrating what still needs to be done. So my suggestion: Would you mind to fix docstrings, check for ndimensionality of the code I suggested and add your suggested changes for passed elements and privacy of the functions? Once it's online I will check it and add what needs to be done. What do you think?

[the-lay]

Hey, I've finally got around to work on tiler. Sorry it takes me so long :(

There were a lot of changes in master branch, so I merged it here and also did some new general changes (latter by mistake, but that's fine since I still want to merge this back to master eventually). I am now looking over your code/PR more thoroughly and want to finish this weekend.

Hence: division by the weight map should be the right solution. what do you think?

I agree! But that shouldn't be default, rather a flag parameter in Merger.merge IMO.

Would you mind to fix docstrings, check for ndimensionality of the code I suggested and add your suggested changes for passed elements and privacy of the functions?

Sure, will focus on that now.

A question from me: What do you mean by minimal overlap, could you explain Tiler.calculate_minimal_overlap method? After all, even overlap of 0 is fine, depends on the use cases.

[renerichter]

Hiho, no worries. I wanted to answer already two weeks ago but somehow didn't manage. Thank you for the merge and comments.

I agree! But that shouldn't be default, rather a flag parameter in Merger.merge IMO. Why would this be better. I believe that people just want to tile an object, calculate whatever on each tile and then merge it back together with a neat overleap. Hence: a flag parameter that is set to active on default should be the way to go, or?

What do you mean by minimal overlap, could you explain Tiler.calculate_minimal_overlap method? After all, even overlap of 0 is fine, depends on the use cases.

Actually, the function Tiler.calculate_minimal_overlap is a remnant of my approach towards the automated tiling calculation and I thought that having such a method for a given "total image size" + " amount of tiles to be split into" would be helpful. Then I decided (like was basically in the toolbox already) that rather having the user decide how big a tile (and the overlap) should be and from there on calculating amount and positions of tiles is more useful and straight forward. Hence: We can delete it.

On another note: CI are active for this repo. Did you test and integrate them properly? Shall we have a further look and or test?

[the-lay]

I have pushed a big pack of commits. I have refactored your code to be a bit more end user friendly and consistent, hope you don't mind. I've added/modified tests, documentation, fixed CI on pull requests, and generally added few new features (sorry for doing that in your fork, but hopefully I can merge it back very soon!).

Now padding is applied through Tiler.calculate_padding(), see the example here. I have also worked on normalization by weight. I think now it's a bit less complex, just a matter of setting flag normalize_by_weights=True in Merger.merge().

Here's the adjusted code to your code examples and figures with the new changes:

Under spoiler

```python import numpy as np from tiler import Tiler, Merger import matplotlib.pyplot as plt # %% Test example for online # parameter data_shape = [1, 256, 256] tile_shape = np.array([1, 64, 64]) overlap = np.array([0, 7, 20]) tiler_mode = 'wrap' windows_supported = ['boxcar', 'triang', 'blackman', 'hamming', 'hann', 'bartlett', 'parzen', 'bohman', 'blackmanharris', 'nuttall', 'barthann', 'overlap-tile'] # image xy = np.ogrid[0:data_shape[-2], 0:data_shape[-1]] xy = np.sqrt(xy[0] * np.transpose(xy[0]) + xy[1] * np.transpose(xy[1])) im = np.cos(4 * xy * np.pi / np.max(xy))[np.newaxis] merged_images = [im, ] # tile me tiler = Tiler(data_shape=data_shape, tile_shape=tile_shape, overlap=overlap, channel_dimension=0) new_shape, padding = tiler.calculate_padding() tiler.recalculate(data_shape=new_shape) im_padded = np.pad(im, padding, mode="reflect") weights_sums = [np.ones(im.shape), ] for mwin in windows_supported: merger = Merger(tiler, window=mwin) for tile_id, tile in tiler(im_padded): processed_tile = tile # lambda x: x merger.add(tile_id, processed_tile) imf = merger.merge(extra_padding=padding) merged_images.append(imf) weights_sums.append(merger._unpad(merger.weights_sum, padding)) # plot images merged_images = np.array(merged_images)[:, 0] fig, ax = plt.subplots(nrows=4, ncols=4, figsize=[14, 14]) plt_titles = ['reference', ] + windows_supported axm = ax.flatten() for m, mim in enumerate(merged_images): ima = axm[m].imshow(mim, interpolation='None') axm[m].set_title(plt_titles[m]) plt.suptitle('Resulting Images') plt.tight_layout() plt.show() # plot images weights_sums = np.array(weights_sums) fig2, ax2 = plt.subplots(nrows=4, ncols=4, figsize=[14, 14]) plt_titles = ['reference', ] + windows_supported axm = ax2.flatten() with np.errstate(divide='ignore', invalid='ignore'): inv_weights = np.nan_to_num(weights_sums[0] / weights_sums) inv_weights /= np.max(inv_weights, axis=(-2, -1), keepdims=True) for m, invm in enumerate(inv_weights): ima = axm[m].imshow(invm[0] ** 0.002, interpolation='None') axm[m].set_title(plt_titles[m]) plt.suptitle('Inverse Weights-maps **0.002 (for display)') plt.tight_layout() plt.show() # differences for m, mim in enumerate(merged_images): print(f"All close for window_func={plt_titles[m]}?\t {np.allclose(mim, im[0])}") ``` ``` All close for window_func=reference? True All close for window_func=boxcar? True All close for window_func=triang? True All close for window_func=blackman? True All close for window_func=hamming? True All close for window_func=hann? True All close for window_func=bartlett? True All close for window_func=parzen? True All close for window_func=bohman? True All close for window_func=blackmanharris? True All close for window_func=nuttall? True All close for window_func=barthann? True All close for window_func=overlap-tile? True ```  

===

Why would this be better. I believe that people just want to tile an object, calculate whatever on each tile and then merge it back together with a neat overleap. Hence: a flag parameter that is set to active on default should be the way to go, or?

After thinking about it some more, I believe you are right, thank you for defending the on point! On the first thought I thought it would not work when window is not specified, but it all works out nicely even then.

Tiler.calculate_minimal_overlap is a remnant of my approach towards the automated tiling calculation and I thought that having such a method for a given "total image size" + " amount of tiles to be split into" would be helpful.

I actually like the idea! I was thinking about creating few Tiler class methods for auto overlap calculation and for auto tile shape calculation.

[renerichter]

Hey, thank you for all your efforts. Now the code looks way nicer. :) Working on my fork is totally fine and thank you for already including the changes in the main-toolbox. Even though it was more or less only a wrapper I liked (for closeness of usage of the class) that applying the padding via im_padded = tiler.pad_outer(im, tiler.pads) was possible in my code suggestion. Now the padding is again applied via im_padded = np.pad(im, padding, mode="reflect"), which is fine as well, but I liked the class-inherent usage a bit more as then the call could be standardized to tiler.pad_outer() which then internally applies tiler.pads to im and stores it (maybe) in another object called tiler._im_padded. Hence the user needs less and less to think about what is happening behind the scenes. What do you think? Again: Thank you for all your good work!

[the-lay]

I see how wrapping np.pad or storing padded image in object can be easier for users, but I think it might be restricting them too much.

For example, in segmentation tasks you can have many "layers" of the same data, e.g. RGB image and class mask. In my workflow, I'd create one Tiler object and apply calculated padding to all images and with different padding modes, essentially making pad_outer probably same to np.pad. Then, obtain tiles of multiple images at once with the same Tiler but passing different data arrays. If there is a im_padded inside, I'd need to create separate Tilers with the same settings.

So I'd prefer to just explicitly return padding in the form that is easy to plug into np.pad and let users handle padded images.

I will merge this all into main branch later today and make a new version soon! Thank you very much for your help and discussions!