snibgo's ImageMagick pages

Gaussian blur by integral images

Fast approximations.

• References
• Sample inputs
• Colorspaces
• Concepts
• Approximation with resize
• Approximation with single 2D deintegration
• Approximation with single 1D deintegration
• Approximation with repeated integral images
• Exact replication with SII
• Approximation with multiple different deintegrations (SII)
• Multi-method script
• Performance
• Future
• Scripts

This page builds on material in Integral images. It is mostly about replicating 1D Gaussian blurs using integral images.

A 2D Gaussian blur is separable: it is numerically equivalent to applying a 1D Gaussian blur horizontally (or vertically), then applying another 1D Gaussian blur in the other direction.

Scripts on this page assume that the version of ImageMagick in %IMDEV% has been built with various process modules (see Process modules), and is HDRI.

References

This page was inspired by:

• Stacked Integral Image Amit Bhatia, Wesley E. Snyder and Griff Bilbro, 2010.
• Efficient and Accurate Gaussian Image Filtering Using Running Sums Elhanan Elboher and Michael Werman, 2011.
• A Survey of Gaussian Convolution Algorithms Pascal Getreuer, 2013.

The first reference is vague on how to calculate the sizes and weights, but a method is given in the other two references.

See also:

• Fast Almost-Gaussian Filtering Peter Kovesi, 2010.

Sample inputs

For convenience, most examples use a small image, with small blurs.

set SRC=toes.png

However, the benefits are seen when we use large images and large blurs.

set LGE_SRC=zp_sus.tiff

%IMG7%magick %LGE_SRC% info: 

zp_sus.tiff TIFF 4924x7378 4924x7378+0+0 16-bit TrueColor sRGB 207.88MiB 0.000u 0:00.000

We will show a small crop from this large image.

set LGE_CROP=-crop 520x440+1593+3455 +repage set LGE_SIG=100 %IMG7%magick ^ %LGE_SRC% ^ %LGE_CROP% ^ gii_lge_crp.png

Colorspaces

set ColSpIn=-colorspace RGB -set colorspace sRGB
set ColSpOut=-set colorspace RGB -colorspace sRGB

set ColSpIn=
set ColSpOut=

Concepts

A Gaussian blur is a convolution with a kernel that is theoretically infinitely large, so every output pixel is influenced by every input pixel. The influence (or "weight") of distant pixels is less than the influence of nearby pixels. In practice, when IM calculates the kernel, it ignores weights that are smaller than a certain threshold, so the kernel is not infinitely large; it is truncated. The threshold depends on the Q-number of the IM build. We demonstrate this with two builds of IM:

%IMG7%magick ^
  xc: ^
  -format "Q%%q\n" ^
  -define morphology:showKernel=1 ^
  -morphology Convolve:0 "Blur:0x0.4;Blur:0x0.4,90" ^
  info: >gii_knl1.lis 2>&1

Kernel #0 "Blur" of size 3x1+1+0 with values from 0 to 0.802398
Forming a output range from 0 to 1 (Normalized)
 0: 0.0988012  0.802398 0.0988012
Kernel #1 "Blur@90" of size 1x3+0+1 with values from 0 to 0.802398
Forming a output range from 0 to 1 (Normalized)
 0: 0.0988012
 1:  0.802398
 2: 0.0988012
Q16

%IM7DEV%magick ^
  xc: ^
  -format "Q%%q\n" ^
  -define morphology:showKernel=1 ^
  -morphology Convolve:0 "Blur:0x0.4;Blur:0x0.4,90" ^
  info: >gii_knl2.lis 2>&1

Kernel #0 "Blur" of size 5x1+2+0 with values from 0 to 0.802305
Forming a output range from 0 to 1 (Normalized)
 0: 5.77217e-05 0.0987898  0.802305 0.0987898 5.77217e-05
Kernel #1 "Blur@90" of size 1x5+0+2 with values from 0 to 0.802305
Forming a output range from 0 to 1 (Normalized)
 0: 5.77217e-05
 1: 0.0987898
 2:  0.802305
 3: 0.0987898
 4: 5.77217e-05
Q32

For the same sigma=0.4, the kernel for Q16 has three elements, and for Q32 it has five elements.

This page assumes that "-blur N" and "-morphology Convolve Blur:0xN;Blur:0xN,90" use the same kernels and make the same results. We show a summary of the kernels:

%IMG7%magick xc: -define morphology:showKernel=1 -blur 0x10 info: 2>&1 |grep Kernel 

Kernel #0 "Blur" of size 79x1+39+0 with values from 0 to 0.0398826
Kernel #1 "Blur@90" of size 1x79+0+39 with values from 0 to 0.0398826

%IMG7%magick xc: -define morphology:showKernel=1 -morphology Convolve:0 Blur:0x10 info: 2>&1 |grep Kernel 

Kernel "Blur" of size 79x1+39+0 with values from 0 to 0.0398826

And we show the images from these two methods of Gaussian blurring:

%IMG7%magick ^ %SRC% ^ %ColSpIn% ^ -blur 0x10 ^ %ColSpOut% ^ gii_q1.png We will compare results to this image.
%IM7DEV%magick ^ %SRC% ^ %ColSpIn% ^ -morphology Convolve "Blur:0x10;Blur:0x10,90" ^ %ColSpOut% ^ gii_q2.png

The images look the same. How similar are they?

%IMG7%magick ^
  gii_q1.png gii_q2.png ^
  -metric RMSE -format %%[distortion] -compare ^
  info: 

1.35784e-05

They are practically identical, demonstrating that "-blur" and "-morphology Convolve Blur" do the same job.

We show a crop of the large image blurred.

set MAGICK_OCL_DEVICE=true %IM7DEV%magick ^ %LGE_SRC% ^ %ColSpIn% ^ -blur 0x100 ^ -define "quantum:format=floating-point" ^ -depth 32 ^ %ColSpOut% ^ +write gii_lge_q1.tiff ^ %LGE_CROP% ^ gii_lge_q1.png set MAGICK_OCL_DEVICE=false

The file gii_lge_q1.tiff will be used for comparisons.

With script knlWidths.bat, using IM Q16, we can tabulate the width (and height) of the kernel for various values of sigma:

The kernel width is an upper limit on the window size of deintegrations. Windows larger than this are certain to be useless.

The operation "-process deintegim window WxH" sets pixels to the mean of the sliding window. Internally, it calculates the value of the sum of the pixel values, and divides this by the sum of the alpha values. The dontdivide or dd option doesn't divide by anything, so sets pixels to the sum of the sliding window.

Approximation with resize

For comparison with other methods, we blur with "-resize". The resize percentage is chosen to give the result that is closest to the result from "-blur 0x10".

%IMG7%magick ^ %SRC% ^ %ColSpIn% ^ -set option:MYSIZE %%[fx:w]x%%[fx:h] ^ -resize 8.5%% ^ -resize "%%[MYSIZE]^!" ^ %ColSpOut% ^ gii_res.png

The result is clearly different to the result from "-blur 0x10". What is the RMSE?

%IMG7%magick ^
  gii_q1.png gii_res.png ^
  -gravity Center -crop 80x80%%+0+0 +repage ^
  -metric RMSE -format %%[distortion] -compare ^
  info: 

0.0356201

The large image:

%IM7DEV%magick ^ %LGE_SRC% ^ %ColSpIn% ^ -set option:MYSIZE %%[fx:w]x%%[fx:h] ^ -resize 0.47%% ^ -resize "%%[MYSIZE]^!" ^ -define "quantum:format=floating-point" ^ -depth 32 ^ %ColSpOut% ^ +write gii_lge_res.tiff ^ %LGE_CROP% ^ gii_lge_res_sm.png

%IM7DEV%magick ^
  gii_lge_q1.tiff ^
  gii_lge_res.tiff ^
  -gravity Center -crop 80x80%%+0+0 +repage ^
  -metric RMSE -format %%[distortion] -compare ^
  info: 

0.0118526

Approximation with single 2D deintegration

"integim" followed by a single 2D "deintegim window WWxHH" sets each output pixel to the mean of a rectangle of input pixels. This is a blur, and we can measure how close it is to a "-blur" operation. The script iiGauss2d.bat finds what square window diameter BestW is closest to a blur of a given sigma. We tabulate results for some values of sigma:

call %PICTBAT%iiGauss2dN gii_2d.htm

An algorithm for this is:

calcR = max (1, floor (sigma * 1.63) );
if (sigma > 100) {
  calcR -= floor (0.5 + 575 * pow ( (sigma-100)/1900, 1.409 ) );
}
calcW = 2 * calcR + 1;

We show the result from a 2D deintegration:

%IM7DEV%magick ^ %SRC% ^ -process integim ^ -process 'deintegim window 33x33' ^ gii_2d.png

Differences can be seen at the edges, but the result is visually acceptable.

%IMG7%magick ^
  gii_q1.png gii_2d.png ^
  -gravity Center -crop 80x80%%+0+0 +repage ^
  -metric RMSE -format %%[distortion] -compare ^
  info: 

0.00701666

Aside from the edges, the RMSE is within my usual no-visible-difference threshold of 1%.

Approximation with single 1D deintegration

For a given sigma, what single "-process integim -process 'deintegim window Nx1'" gives the closest result to a 1D Gaussian blur?

For a sigma, we can compare a blurred impulse image to images made by integim and deintegim at different window widths, to find the window width that has the smallest RMSE difference. The script iiGaussOne.bat does this by searching odd-number window widths. (A binary chop would be faster, but the search is needed only once.)

call %PICTBAT%iiGaussOne 1

echo sigma=%igoSigma% MaxW=%igoMaxW% BestW=%igoBestW% 

sigma=1 MaxW=9 BestW=3 

We repeat this for different values of sigma:

call %PICTBAT%iiGaussOneN gii_g1n.htm

As we might expect, the best deintegration width is nearly proportional to the required sigma: roughly BestW = 2 * floor (1.53 * sigma) + 1. But how good is it?

%IM7DEV%magick ^ %SRC% ^ -process integim ^ -process 'deintegim window 31x1' ^ -process integim ^ -process 'deintegim window 1x31' ^ gii_g1bd.png

Compared to "-blur 0x10", differences are visible at the edges, especially towards the right of the bottom edge. But the result is visually acceptable.

%IMG7%magick ^
  gii_q1.png gii_g1bd.png ^
  -gravity Center -crop 80x80%%+0+0 +repage ^
  -metric RMSE -format %%[distortion] -compare ^
  info: 

0.00667103

Aside from the edges, the RMSE is within my normal acceptable limit of 1%. The RMSE is slightly better than for the single 2D deintegration, but the processing has doubled.

[Simpler method: ??]

Approximation with repeated integral images

The previous section showed how a single "-process integim -process 'deintegim window Nx1'" closely approximates a 1D Gaussian blur.

But we can repeat the integral images processes (integim and deintegim), for example three times, using the same parameters each time. For a given sigma, what is the best window width?

The script iiGaussMult.bat searches odd-number window widths up to MaxW, where MaxW is the width of the blur kernel.

We tabulate results for 3 iterations:

call %PICTBAT%iiGaussMultN gii_gmultn.htm

Values of BestW are smaller than in the previous section. Again, BestW is roughly proportional to sigma: BestW = 2*floor (0.95 * sigma) + 1.

This needs smaller deintegration windows than the previous section. This is no advantage, as processing time is independent of window size.

Using a window width of 19, we make a 2D blur:

%IM7DEV%magick ^ %SRC% ^ -process integim ^ -process 'deintegim window 19x1' ^ -process integim ^ -process 'deintegim window 19x1' ^ -process integim ^ -process 'deintegim window 19x1' ^ -process integim ^ -process 'deintegim window 1x19' ^ -process integim ^ -process 'deintegim window 1x19' ^ -process integim ^ -process 'deintegim window 1x19' ^ gii_g1bdm.png

%IMG7%magick ^
  gii_q1.png gii_g1bdm.png ^
  -gravity Center -crop 80x80%%+0+0 +repage ^
  -metric RMSE -format %%[distortion] -compare ^
  info: 

0.00189328

The result is excellent, a five-fold improvement in RMSE, but needing three times the amount of processing.

What happens if we interleave the processes?

%IM7DEV%magick ^ %SRC% ^ -process integim ^ -process 'deintegim window 19x1' ^ -process integim ^ -process 'deintegim window 1x19' ^ -process integim ^ -process 'deintegim window 19x1' ^ -process integim ^ -process 'deintegim window 1x19' ^ -process integim ^ -process 'deintegim window 19x1' ^ -process integim ^ -process 'deintegim window 1x19' ^ gii_g1bdmi.png

%IMG7%magick ^
  gii_q1.png gii_g1bdmi.png ^
  -gravity Center -crop 80x80%%+0+0 +repage ^
  -metric RMSE -format %%[distortion] -compare ^
  info: 

0.00189328

The result is the same, so interleaving has no advantage.

Is there anything magical about three? How about repeating ten times? For ten iterations, we can calculate that BestW=11.

%IM7DEV%magick ^ %SRC% ^ -process integim ^ -process 'deintegim window 11x1' ^ -process integim ^ -process 'deintegim window 11x1' ^ -process integim ^ -process 'deintegim window 11x1' ^ -process integim ^ -process 'deintegim window 11x1' ^ -process integim ^ -process 'deintegim window 11x1' ^ -process integim ^ -process 'deintegim window 11x1' ^ -process integim ^ -process 'deintegim window 11x1' ^ -process integim ^ -process 'deintegim window 11x1' ^ -process integim ^ -process 'deintegim window 11x1' ^ -process integim ^ -process 'deintegim window 11x1' ^ -process integim ^ -process 'deintegim window 1x11' ^ -process integim ^ -process 'deintegim window 1x11' ^ -process integim ^ -process 'deintegim window 1x11' ^ -process integim ^ -process 'deintegim window 1x11' ^ -process integim ^ -process 'deintegim window 1x11' ^ -process integim ^ -process 'deintegim window 1x11' ^ -process integim ^ -process 'deintegim window 1x11' ^ -process integim ^ -process 'deintegim window 1x11' ^ -process integim ^ -process 'deintegim window 1x11' ^ -process integim ^ -process 'deintegim window 1x11' ^ gii_g1bdmi2.png

%IMG7%magick ^
  gii_q1.png gii_g1bdmi2.png ^
  -gravity Center -crop 80x80%%+0+0 +repage ^
  -metric RMSE -format %%[distortion] -compare ^
  info: 

0.00046896

The result from ten iterations is four times better than with three iterations, but increases processing by a factor of 3.33.

The Getreuer paper, citing Wells, sugests calculating r (the radius of the box filter) from K (usually 3, 4 or 5) and sigma (σ):

σ2 = K/12 * ((2r+1)2 - 1)

r = floor (0.5 * sqrt(12/K * σ2 + 1))

Defining the window diameter D = 2r+1:

σ2 = K/12 * (D2 - 1)

D = 1 + 2 * floor (0.5 * sqrt(12/K * σ2 + 1))

From that, we can tabulate values of D from K and sigma, according to Wells:

call %PICTBAT%iiWells.bat 

The column for K=3, where roughly D = 2*sigma + 1, can be compared to my finding above, where sigma=1000 has BestW=1899.

Note that for large sigma, we can simplify the Wells formula to:

D = floor (sqrt(12/K) * σ)+0.5)

When K=3, this is:

D = floor (2 * σ) + 0.5

%IM7DEV%magick ^ %SRC% ^ -process integim ^ -process 'deintegim window 21x1' ^ -process integim ^ -process 'deintegim window 21x1' ^ -process integim ^ -process 'deintegim window 21x1' ^ -process integim ^ -process 'deintegim window 1x21' ^ -process integim ^ -process 'deintegim window 1x21' ^ -process integim ^ -process 'deintegim window 1x21' ^ gii_g1bdm2.png

%IMG7%magick ^
  gii_q1.png gii_g1bdm2.png ^
  -gravity Center -crop 80x80%%+0+0 +repage ^
  -metric RMSE -format %%[distortion] -compare ^
  info: 

0.00426322

This result with width=21 compares to my result from width=19 above. The result with 21, though still acceptable, is worse, so I will stick to my formula.

This result is good, but expensive, requiring three repetitions of integim and deintegim. We can get a better RMSE at a lower cost.

Exact replication with SII

Before considering a lower cost approximation, we will consider a perfect replication of -blur using stacked integral images.

A 1D horizontal Gaussian blur can be approximated by calculating the integral image, and deintegrating that with a number of windows at different widths (and constant height of one), and summing those deintegrals with suitable weights. This is the Stacked Integral Images, SII, method.

To explain the process, consider a very simple example: a convolution with the kernel [0.1,0.8,0.1]. This kernel sums to one. Each output pixel value will be set to 0.8 times the corresponding input, plus 0.1 times the value to the left, plus 0.1 times the value to the right.

outx = 0.1*inx-1 + 0.8*inx + 0.1*inx+1 

We can draw a graph of the kernel:

Put another way, each output is 0.7 times the corresponding input, plus 0.1 times the sum (the sum, not the mean) of three input pixels. "0.7" is calculated from the weight of the kernel element minus the element on either side.

outx = 0.7*inx + 0.1*( inx-1 + inx + inx+1 ) 

So instead of convolving with the variable kernel [0.1,0.8,0.1], we convolve with two kernels [0.7] and [0.1,0.1,0.1], and add the results.

Even better, we can convolve with two kernels [1] and [1,1,1], making result images R1 and R2, then calculate R1*0.7+R2*0.1. The convolutions can be performed by making a single integral image which is deintegrated multiple times, and we stack the results, hence the term Stacked Integral Images.

Also note that convolving with kernel [1] has no effect, so we don't need to do that convolution.

So we can exactly replicate a 1D Gaussian blur of size 2*N+1 by creating an integral image and deintegrating it N-1 times. A Gaussian blur with a kernel of 101 elements would need 49 deintegrations.

We show a more realistic example, blurring an ordinary photo.

%IMG7%magick ^ %SRC% ^ -morphology Convolve Blur:0x1 ^ gii_blr1.png

%IMG7%magick ^
  xc: ^
  -define morphology:showKernel=1 ^
  -morphology Convolve:0 Blur:0x1 ^
  NULL: 

Kernel "Blur" of size 9x1+4+0 with values from 0 to 0.384572
Forming a output range from 0 to 1 (Normalized)
 0: 0.000215825 0.00579007  0.059897  0.241811  0.384572  0.241811  0.059897 0.00579007 0.000215825

These factors sum to one.

We can graph the kernel, by convolving an impulse image (one white pixel surrounded by black pixels).

%IMG7%magick ^ -size 100x1 xc:Black ^ -size 1x1 xc:White ^ -size 100x1 xc:Black ^ +append +repage ^ -morphology Convolve Blur:0x1 ^ -bordercolor Black -border 1 ^ -trim +repage ^ -scale "600x1^!" ^ gii_imp1.png call %PICTBAT%graphLineCol ^ gii_imp1.png . . 0 . . . ^ "-compose Copy -bordercolor #888 -border 5"

The area under the graph is one. The kernel is symmetrical, with width=9, so we can replicate the effect with five stacked integral images. To calculate the five weights from the 9x1 kernel, we use the five final values in the kernel, and calculate each weight as a kernel value minus the next kernel value.

If we used the dontdivide option the edges of the image would look terrible. This is because the deintegration takes the sum of the window, but the windows are smaller at the edges. This would mess with the "-poly" calculation; fewer values would be used, so the total would be smaller than it would be nearer the centre. Instead, we use the adjedges (ae) option which multiples by requested_window_size / actual_window_size.

rem Diff blur set W1=0.142761 set W2=0.181914 set W3=0.05410693 set W4=0.005574245 set W5=0.000215825 %IM7DEV%magick ^ %SRC% +write mpr:SRC ^ -process integim ^ mpr:SRC ^ ( -clone 0 -process 'deintegim window 3x1 ae' ) ^ ( -clone 0 -process 'deintegim window 5x1 ae' ) ^ ( -clone 0 -process 'deintegim window 7x1 ae' ) ^ ( -clone 0 -process 'deintegim window 9x1 ae' ) ^ -delete 0 ^ -poly "%W1%,1 %W2%,1 %W3%,1 %W4%,1 %W5%,1" ^ gii_blr2.png

How close is this to the "-morphology Convolve Blur:0x1" operation, comparing just the central 80%?

%IMG7%magick ^
  gii_blr1.png gii_blr2.png ^
  -gravity Center -crop 80x80%%+0+0 +repage ^
  -metric RMSE -format %%[distortion] -compare ^
  info: 

1.38446e-06

The RMSE is practically zero, so the results are practically identical.

As the general method is called stacked integral images, we might call this exact method a full stack.

Beware that this method can eat memory. If the kernel radius is 500, this method will need 501 copies of the image in memory simultaneously for the "-poly" operation. A modified version could be written to add weighted images as they are made. This would take longer, but wouldn't eat memory.

The script iiGaussExact.bat implements this method.

call %PICTBAT%iiGaussExact %SRC% gii_exact.png 10

How close is this to "-blur 0x10", apart from the image edges?

%IMG7%magick ^
  gii_q1.png gii_exact.png ^
  -gravity Center -crop 80x80%%+0+0 +repage ^
  -metric RMSE -format %%[distortion] -compare ^
  info: 

0.00015076

They are practically identical.

Approximation with multiple different deintegrations (SII)

[[When we have a stack of K deintegrations, there are 2*K-1 unknowns: the width of each kernel, and the weight of each deintegration, but minus one because we know the area under the graph should be one, ie:

i=k
sum (widthi * weighti) = 1
i=1

When K=1, there is only one deintegration in each direction, and there is one unknown. We have covered that case in Approximation with single 1D deintegration above. ]]

The section discusses two methods for calculating a stack of different deintegrations to approximate a 1D Gaussian blur.

The outline algorithm for K slices, numbered from 0 at the top to K-1 at the bottom, designating the vertical axis as V, is:

for (i=0; i < K, i++) {
  Calculate the Vtop of this slice.
  Decide the required weight of this slice.
  Calculate the Vbottom of this slice.
  Find the radius of this slice.
}
Calculate the total area of the slices.
Adjust weights to make the total area 1.0.

In each slice:

Vtop - Vbottom = weight

For the top slice, V_top is somewhere near the maximum value of the Gaussian curve. For the bottom slice, V_bottom is zero.

How do we "Decide the required weight of this slice"? We might make the weights for all slices equal. By experimentation, we find that reducing the required weight of the bottom slice, and increasing the weights of the other slices to compensate, is a slight improvement.

The two methods differ in how they "Find the radius of this slice". Ideally, the radius is such that the RMSE between the block and the slice of the curve (each being a Nx1 image) is minimised. Each value in the block is the same as the maximum value in the slice of the curve.

Mid-edge method

The radius of each slice is calculated such that the value of the Gaussian curve at the slice edge, V, is (V_top+V_bottom)/2. The curve intersects each slice at the mid-height, like these invented examples:

In each horizontal slice, the coloured block has a three-sided area that is outside the curve so the block over-estimates, and the curve has a three-sided area that is outside the coloured block so the block under-estimates. These areas are approximately equal, so they roughly cancel each other. Where the curve is actually a straight line, the three-sided areas will be exactly equal.

The script iiGaussKcalc.bat does the calculations:

call %PICTBAT%iiGaussKcalc 10 3 gii_kcalc_

set gii_kcalc_ 

gii_kcalc_Diam[0]=13
gii_kcalc_Diam[1]=27
gii_kcalc_Diam[2]=41
gii_kcalc_K=3
gii_kcalc_Sigma=10
gii_kcalc_Weight[0]=0.0148514851485149
gii_kcalc_Weight[1]=0.0148514851485149
gii_kcalc_Weight[2]=0.00990099009900987

Now we can transform an image using the calculations: call %PICTBAT%iiGaussK %SRC% gii_kcalc1.png gii_kcalc_

The result is visually good. How close is this to the "-blur 0x10" result?

%IMG7%magick ^
  gii_q1.png gii_kcalc1.png ^
  -gravity Center -crop 80x80%%+0+0 +repage ^
  -metric RMSE -format %%[distortion] -compare ^
  info: 

0.00338361

Slice area method

The radius of each slice is calculated such that when viewed as a graph, the area of the block (diameter*weight) is equal to the area of the Gaussian curve between V_top and V_bottom. This is not totally trivial.

The area of one side of a slice of the curve is calculated easily: cap the values at Vtop, and subtract Vbottom. (But for the top slice, we don't cap values.) Dividing that by the required weight gives the radius, a floating-point number, and we can calculate diameter=radius*2+1.

But the radius must be an integer, so the areas will probably not be exactly equal. We adjust the weight of the block to make the areas equal, and this adjusts the value of V_top. But now a different part of the curve is relevant, so we need to re-calculate the area of the curve, and so on. We iterate until the relative difference between the areas is below a desired threshold (the error limit).

The script iiGaussKcalc_4.bat calculates the deintegration windows and weights from the required sigma and the number of slices and the error limit, and the script iiGaussK.bat blurs an image from the calculation.

The number of slices is from 1 upwards. If only one slice is needed, the method Approximation with single 1D deintegration above, which explores different weights and hence different diameters, will give better results.

An error limit of 1.0 won't repeat the calculation to adjust the areas. An error limit of 0.0 will repeat until the areas are equal, to 15 decimal places. With rounding errors this may never happen, so there is a hard limit of 100 iterations. A limit of around 1e-8 is generally low enough for high-precision results. If we are doing the calculation just once to blur many images, we don't care how long the calculation takes, so an error limit of 0 may be appropriate.

Do the calculations:

call %PICTBAT%iiGaussKcalc_4 10 3 gii_gk_

The calculations are recorded as environment variables, and also in a file. In this case, the file is named gii_gk__vars.lis.

gii_gk_Diam[0]=13
gii_gk_Diam[1]=27
gii_gk_Diam[2]=45
gii_gk_K=3
gii_gk_Sigma=10
gii_gk_Weight[0]=0.014503254618265
gii_gk_Weight[1]=0.016733613815784
gii_gk_Weight[2]=0.00800631025408104

Now we can transform an image using the calculations: call %PICTBAT%iiGaussK %SRC% gii_gk1.png gii_gk_

The result is visually good. How close is this to the "-blur 0x10" result?

%IMG7%magick ^
  gii_q1.png gii_gk1.png ^
  -gravity Center -crop 80x80%%+0+0 +repage ^
  -metric RMSE -format %%[distortion] -compare ^
  info: 

0.00203608

The result is substantially better than the mid-edge method.

The script iiGaussKgraph.bat draws a graph of the situation. It shows the Gaussian curve as a white line, and each deintegration as a coloured block. The whole is scaled vertically, because large sigmas make very flat curves.

call %PICTBAT%iiGaussKgraph ^ gii_gk_ gii_gk1_gr.png

Another script, iiGaussChkRmse.bat, re-calculates slices from given weights by finding the block radius that minimises RMSE from the slice of the curve. (However, it does not adjust weights for the integer radius problem.)

call %PICTBAT%iiGaussChkRmse gii_gk_ gii_recalc_

set gii_recalc_ 

gii_recalc_Diam[0]=15
gii_recalc_Diam[1]=29
gii_recalc_Diam[2]=45
gii_recalc_K=3
gii_recalc_Sigma=10
gii_recalc_Weight[0]=0.014503254618265
gii_recalc_Weight[1]=0.016733613815784
gii_recalc_Weight[2]=0.00800631025408104

The graph and the blurred image look like this:

call %PICTBAT%iiGaussKgraph ^ gii_recalc_ gii_recalc_gr.png
call %PICTBAT%iiGaussK %SRC% gii_recalc.png gii_recalc_

How close is that result to -blur 0x10?

%IMG7%magick ^
  gii_q1.png gii_recalc.png ^
  -gravity Center -crop 80x80%%+0+0 +repage ^
  -metric RMSE -format %%[distortion] -compare ^
  info: 

0.0674937

The result is not good.

We show a range of parameters.

call %PICTBAT%iiGaussKsamps toes.png gii_samps

Parameter Sigma	Parameter K	Image	RMSE
10	1		0.0149959
10	2		0.00425315
10	3		0.00203608
10	4		0.0012225
10	5		0.00120234
10	10		0.000305671

Future: The script relies on symmetry in the 1D kernel. It could be modified to remove that reliance -- it would process entire slices instead of half of each slice. The script starts by blurring an impulse image to obtain the kernel it is trying to approximate. It could instead use a kernel as input, which would allow approximation of any 1D kernel. Extending the method to a general non-seperable 2D kernel is more complex. Each slice of that kernel is an arbitrary 2D shape, and we would need to find the rectangle of the same area that best represents that shape. Perhaps this is best done with moments.

Multi-method script

For general-purpose, we need a script iiGaussOpn.bat that creates a suitable environment variable with the required integim etc. Takes parameters: sigma, method, etc.

Performance

The motivation behind this page is to find fast methods to approximate Gaussian blurs of large images with large sigmas.

Future

Some algorithms shown here make an integral image that is deintegrated only in the horizontal direction, then another integral image that is deintegrated only in the vertical direction. The integim module is general-purpose, and can be deintegrated in either direction, or both. Specialist versions could be written for deintegration in only one direction. These would be slightly faster, and more easily parallelised.

Scripts

For convenience, .bat scripts are also available in a single zip file. See Zipped BAT files.

knlWidths.bat

@echo off

echo ^<table^>

echo ^<tr^>^<th^>sigma^</th^>^<th^>width^</th^>^</tr^>

for %%S in (0.6, 1, 2.5, 5, 10, 50, 100, 200, 500, 1000, 2000) do (

  for /F "usebackq tokens=1,5 delims=x+ " %%A in (`%IMG7%magick xc: ^
    -define "morphology:showKernel=1" ^
    -morphology Convolve:0 Blur:0x%%S ^
    NULL: 2^>^&1`) do (
    if "%%A"=="Kernel" echo ^<tr^>^<th^>%%S^</th^>^<td^>%%B^</td^>^</tr^>
  )
)

echo ^</table^>

iiGauss2d.bat

rem What single integral image is closest to a 2D Gaussian blur?

rem %1 Sigma

setlocal enabledelayedexpansion

set sig=%1
if "%sig%"=="." set sig=
if "%sig%"=="" set sig=10

for /F "usebackq" %%L in (`%IMG7%magick ^
  -size %%[fx:%sig%*10]x1 xc:Black ^
  -size 1x1 xc:White ^
  ^( -clone 0 ^) ^
  +append +repage ^
  -define "quantum:format=floating-point" -depth 32 ^
  -morphology Convolve Blur:0x%sig% ^
  -bordercolor Black -border 1 ^
  -trim +repage ^
  -format "MaxW=%%w\n" ^
  info:`) do set %%L

set /A nBord=(%MaxW%-1)/2

echo %0: MaxW=%MaxW% nBord=%nBord%  Blurring...

set SAVE_OCL=MAGICK_OCL_DEVICE
set MAGICK_OCL_DEVICE=true

%IMG7%magick ^
  -size 1x1 xc:White ^
  -bordercolor Black -border %nBord% ^
  -define "quantum:format=floating-point" -depth 32 ^
  +write igo_imp.miff ^
  -blur 0x%sig% ^
  igo_blr.miff

set MAGICK_OCL_DEVICE=%SAVE_OCL%

echo %0: Blurring ended.

%IM7DEV%magick ^
  igo_imp.miff ^
  -process 'integim v' ^
  -define "quantum:format=floating-point" -depth 32 ^
  igo_integ.miff

set MinDiff=1.0
set BestW=
set Bust=0

:: The process is slow, so start at somewhere plausible.
for /F "usebackq" %%L in (`%IMG7%magick ^
  xc: ^
  -format "StartMax=%%[fx:floor(%sig%*1.4)*2+1]" ^
  info:`) do set %%L


echo off

for /L %%N in (%StartMax%,2,%MaxW%) do (
  if !Bust!==0 (
    for /F "usebackq" %%L in (`%IM7DEV%magick ^
      igo_integ.miff ^
      -process 'deintegim window %%Nx%%N' ^
      igo_blr.miff ^
      -metric RMSE ^
      -precision 15 ^
      -format "Diff=%%[distortion]\n" ^
      -compare ^
      info:`) do set %%L

    for /F "usebackq" %%L in (`%IM7DEV%magick ^
      xc: ^
      -format "Better=%%[fx:!Diff!<!MinDiff!?1:0]\n" ^
      info:`) do set %%L

    echo %0: N=%%N Diff=!Diff! Better=!Better!

    if !Better!==1 (
      set MinDiff=!Diff!
      set BestW=%%N
    ) else (
      set Bust=1
    )
  )

)

rem %IMG7%magick igo_integ.miff -verbose info:

echo %0: sig=%sig% MaxW=%MaxW% MinDiff=%MinDiff% BestW=%BestW%

endlocal && set igoSigma=%sig%&&set igoMaxW=%MaxW%&&set igoBestW=%BestW%

iiGauss2dN.bat

rem %1 output HTM file

echo off

(
  echo ^<table^>
  echo ^<tr^>^<th^>sigma^</th^>^<td^>BestW^</td^>^</tr^>
) >%1

for %%N in (0.6,1,2.5,5,10,20,50,100,200,500,1000,2000) do (
  call %PICTBAT%iiGauss2d %%N
  echo ^<tr^>^<th^>%%N^</th^>^<td^>!igoBestW!^</td^>^</tr^> >>%1
)

echo ^</table^> >>%1


type %1

iiGaussOne.bat

rem What single integral image is closest to a 1D Gaussian blur?

rem %1 Sigma

setlocal enabledelayedexpansion

set sig=%1
if "%sig%"=="." set sig=
if "%sig%"=="" set sig=10

for /F "usebackq" %%L in (`%IMG7%magick ^
  -size %%[fx:%sig%*10]x1 xc:Black ^
  -size 1x1 xc:White ^
  ^( -clone 0 ^) ^
  +append +repage ^
  -define "quantum:format=floating-point" -depth 32 ^
  +write igo_imp.miff ^
  -morphology Convolve Blur:0x%sig% ^
  +write igo_blr.miff ^
  -bordercolor Black -border 1 ^
  -trim +repage ^
  -format "MaxW=%%w\n" ^
  info:`) do set %%L

echo MaxW=%MaxW%

%IM7DEV%magick ^
  igo_imp.miff ^
  -process 'integim v' ^
  -define "quantum:format=floating-point" -depth 32 ^
  igo_integ.miff

set MinDiff=1.0
set BestW=
set Bust=0

:: The process is slow, so start at somewhere plausible.
for /F "usebackq" %%L in (`%IMG7%magick ^
  xc: ^
  -format "StartMax=%%[fx:2*floor(%sig%*1.25)+1]" ^
  info:`) do set %%L

echo off

for /L %%N in (%StartMax%,2,%MaxW%) do (
  if !Bust!==0 (
    for /F "usebackq" %%L in (`%IM7DEV%magick ^
      igo_integ.miff ^
      -process 'deintegim window %%Nx1' ^
      igo_blr.miff ^
      -metric RMSE ^
      -format "Diff=%%[distortion]\n" ^
      -compare ^
      info:`) do set %%L

    for /F "usebackq" %%L in (`%IM7DEV%magick ^
      xc: ^
      -format "Better=%%[fx:!Diff!<!MinDiff!?1:0]\n" ^
      info:`) do set %%L

    echo N=%%N Diff=!Diff! Better=!Better!

    if !Better!==1 (
      set MinDiff=!Diff!
      set BestW=%%N
    ) else (
      set Bust=1
    )
  )
)

rem %IMG7%magick igo_integ.miff -verbose info:

echo sig=%sig% MaxW=%MaxW% MinDiff=%MinDiff% BestW=%BestW%

endlocal && set igoSigma=%sig%&&set igoMaxW=%MaxW%&&set igoBestW=%BestW%

iiGaussOneN.bat

rem %1 output HTM file

echo off

(
  echo ^<table^>
  echo ^<tr^>^<th^>sigma^</th^>^<td^>BestW^</td^>^</tr^>
) >%1

for %%N in (0.6,1,2.5,5,10,20,50,100,200,500,1000,2000) do (
  call %PICTBAT%iiGaussOne %%N
  echo ^<tr^>^<th^>%%N^</th^>^<td^>!igoBestW!^</td^>^</tr^> >>%1
)

echo ^</table^> >>%1


type %1

iiGaussMult.bat

rem What multiple iterations of integral image is closest to a 1D Gaussian blur?

rem %1 Sigma
rem %2 numIter

setlocal enabledelayedexpansion

set sig=%1
if "%sig%"=="." set sig=
if "%sig%"=="" set sig=10

set NumIter=%2
if "%%"=="." set NumIter=
if "%%"=="" set NumIter=10

for /F "usebackq" %%L in (`%IMG7%magick ^
  -size %%[fx:%sig%*10]x1 xc:Black ^
  -size 1x1 xc:White ^
  -size %%[fx:%sig%*10]x1 xc:Black ^
  +append +repage ^
  -define "quantum:format=floating-point" -depth 32 ^
  +write igo_imp.miff ^
  -morphology Convolve Blur:0x%sig% ^
  +write igo_blr.miff ^
  -bordercolor Black -border 1 ^
  -trim +repage ^
  -format "MaxW=%%w\n" ^
  info:`) do set %%L

echo MaxW=%MaxW%

%IM7DEV%magick ^
  igo_imp.miff ^
  -process 'integim v' ^
  -define "quantum:format=floating-point" -depth 32 ^
  igo_integ.miff


set MinDiff=1.0
set BestW=
set Bust=0

:: The process is slow, so start at somewhere plausible.
for /F "usebackq" %%L in (`%IMG7%magick ^
  xc: ^
  -format "StartMax=%%[fx:2*floor(0.8*%sig%)+1]" ^
  info:`) do set %%L

echo off

for /L %%N in (%StartMax%,2,%MaxW%) do if !Bust!==0 (

  set PROC=
  for /L %%I in (1,1,%NumIter%) do set PROC=!PROC! -process integim -process 'deintegim window %%Nx1'

  for /F "usebackq" %%L in (`%IM7DEV%magick ^
    igo_imp.miff ^
    !PROC! ^
    igo_blr.miff ^
    -metric RMSE ^
    -format "Diff=%%[distortion]\n" ^
    -compare ^
    info:`) do set %%L

  for /F "usebackq" %%L in (`%IM7DEV%magick ^
    xc: ^
    -format "Better=%%[fx:!Diff!<!MinDiff!?1:0]\n" ^
    info:`) do set %%L

  echo N=%%N Diff=!Diff! Better=!Better!

  if !Better!==1 (
    set MinDiff=!Diff!
    set BestW=%%N
  ) else (
    set Bust=1
  )

)

rem %IMG7%magick igo_integ.miff -verbose info:

echo sig=%sig% MaxW=%MaxW% MinDiff=%MinDiff% BestW=%BestW%

endlocal && set igoSigma=%sig%&&set igoMaxW=%MaxW%&&set igoBestW=%BestW%

iiGaussMultN.bat

rem %1 output HTM file

echo off

(
  echo ^<table^>
  echo ^<tr^>^<th^>sigma^</th^>^<td^>BestW^</td^>^</tr^>
) >%1

for %%N in (0.6,1,2.5,5,10,20,50,100,200,500,1000,2000) do (
  call %PICTBAT%iiGaussMult %%N 3
  echo ^<tr^>^<th^>%%N^</th^>^<td^>!igoBestW!^</td^>^</tr^> >>%1
)

echo ^</table^> >>%1


type %1

iiGaussExact.bat

rem Replicate a Gaussian blur exactly by integral images,
rem by "fully stacked integral images".
rem
rem %1 input image
rem %2 output image
rem %3 sigma
@rem
@rem Beware: for large values of sigma, this eats memory.
@rem A version of this could be written that adds weighted images as they are made,
@rem instead of keeping each in memory until the "-poly".
@rem

@rem TODO: make environment variables and file.


setlocal enabledelayedexpansion

set TMPDIR=\temp

set fDeinteg1=%TMPDIR%\iige1.txt
set fDeinteg2=%TMPDIR%\iige2.txt
set fPoly=%TMPDIR%\iipoly.txt
set fScr=%TMPDIR%\iige.scr

del %fDeinteg1% 2>nul
del %fDeinteg2% 2>nul
del %fPoly% 2>nul
del %fScr% 2>nul

set INFILE=%1
set OUTFILE=%2

set sigma=%3
if "%sigma%"=="." set sigma=
if "%sigma%"=="" set sigma=10


echo off

echo %0: Calculating kernel

set CNT=0

for /F "usebackq tokens=1,*" %%A in (`%IMG7%magick ^
  xc: ^
  -precision 15 ^
  -define "morphology:showKernel=1" ^
  -morphology Convolve:0 Blur:0x%sigma% ^
  info: 2^>^&1 ^|grep "0:" ^|head -n 1 ^| tr "  " \r\n `) do (
  if not "%%A"=="0:" (
    set KNL_V[!CNT!]=%%A
    set /A CNT+=1
  )
)

rem %CNT% should be an odd number

set KNL_V[


::------

echo %0: CNT=%CNT% Writing subfiles

set KNL_V[!CNT!]=0

set /A nLast=%CNT%-1
set /A nMid=%CNT%/2

echo %0: nMid=%nMid% nLast=%nLast%

rem set KNL_V

for /L %%N in (%nMid%,1,%nLast%) do (
  set /A "Diam=(%%N-%nMid)*2+1"
  rem echo %%N !KNL_V[%%N]!

  call :knlNp1 %%N

  for /F "usebackq" %%L in (`%IMG7%magick ^
    xc: ^
    -precision 15 ^
    -format "Weight=%%[fx:!KNL_V[%%N]!-!vp1!]" ^
    info:`) do set %%L

  echo !Weight!,1,>>%fPoly%
  if !Diam! GTR 1 (
    echo ^( -clone 0 -process 'deintegim window !Diam!x1 adjedges' ^) >>%fDeinteg1%
    echo ^( -clone 0 -process 'deintegim window 1x!Diam! adjedges' ^) >>%fDeinteg2%
  )
)

echo %0: Writing %fScr%
(
  echo -process integim
  echo mpr:SRC
  type %fDeinteg1%
  echo -delete 0
  echo -poly @"%fPoly%"
  echo -process integim
  echo mpr:SRC
  type %fDeinteg2%
  echo -delete 0
  echo -poly @"%fPoly%"
  echo -write %OUTFILE%
  echo -exit
) >%fScr%

echo %0: Executing IM

%IM7DEV%magick %INFILE% +write mpr:SRC -script %fScr%

exit /B 0

::-----------------------------------------------------
:: Subroutine

:knlNp1
set /A np1=%1+1
set vp1=!KNL_V[%np1%]!

exit /B 0

iiGauss.bat

rem Approximate Gaussian blur by Stacked Integral Images (SII).
rem Reference integim.htm#apprxconv
rem Following is based on gaussIIeff.pdf

rem %1 input image
rem %2 output image
rem %3 sigma [default: 10]
rem %4 K (3, 4 or 5) [default: 4]

setlocal enabledelayedexpansion

echo off

set INFILE=%1
set OUTFILE=%2

set sigma=%3
if "%sigma%"=="." set sigma=
if "%sigma%"=="" set sigma=10

set sig0=(100/PI)

set K=%4
if "%K%"=="." set K=
if "%K%"=="" set K=4

if %K%==3 (
  set PIs=23,46,76
  set WIs=0.9495,0.5502,0.1618,0
) else if %K%==4 (
  set PIs=19,37,56,82
  set WIs=0.9649,0.6700,0.3376,0.0976,0
) else if %K%==5 (
  set PIs=16,30,44,61,85
  set WIs=0.9738,0.7596,0.5031,0.2534,0.0739,0
) else (
  echo %0: invalid K [%K%]
  exit /B 1
)

set I=1
for %%N in (%PIs%) do (
  set Pi[!I!]=%%N
  set /A I+=1
)

set I=1
for %%N in (%WIs%) do (
  set Wi[!I!]=%%N
  set /A I+=1
)

rem When sigma==sig0,
rem P[i] = Pi[i] roughly,
rem but W[i] = Wi[i]/2 roughly. This seems wrong.
rem Why are weights adjusted at all? To make the area under the curve sum to one?
rem But the 1D case works by adding 3 (or 4 or 5) weighted averages, so the weights should sum to one.

for /L %%I in (1,1,%K%) do (

  echo     -format "P[%%I]=%%[fx:floor(%sigma%/%sig0%*!Pi[%%I]!)]"

  for /F "usebackq" %%L in (`%IMG7%magick xc: ^
    -precision 16 ^
    -format "P[%%I]=%%[fx:floor(%sigma%/%sig0%*!Pi[%%I]!)]" ^
    info:`) do set %%L

  call :vplus1 %%I

  for /F "usebackq" %%L in (`%IMG7%magick xc: ^
    -precision 16 ^
    -format "W[%%I]=%%[fx:!Pi[%%I]!*!Wi[%%I]! / (2*!P[%%I]!+1)]" ^
    info:`) do set %%L

  for /F "usebackq" %%L in (`%IMG7%magick xc: ^
    -precision 16 ^
    -format "Wadj[%%I]=%%[fx:!Pi[%%I]!*(!Wi[%%I]!-!Wp1!) / (2*!P[%%I]!+1)]" ^
    info:`) do set %%L

  for /F "usebackq" %%L in (`%IMG7%magick xc: ^
    -precision 16 ^
    -format "W[%%I]=%%[fx:!Wi[%%I]!-!Wp1!]" ^
    info:`) do set %%L

)

set Pi[
set P[

set Wi[
set W[
set Wadj[

if %K%==3 (
  set PROCS1=
  set PROCS2=
  set sPOLY=-poly "%W[1]%,1 %W[2]%,1 %W[3]%,1"

) else if %K%==4 (
  set PROCS1=^( -clone 0 -process 'deintegim window %P[4]%x1' ^)
  set PROCS2=^( -clone 0 -process 'deintegim window 1x%P[4]%' ^)
  set sPOLY=-poly "%W[1]%,1 %W[2]%,1 %W[3]%,1 %W[4]%,1"
) else (
  set PROCS1=^( -clone 0 -process 'deintegim window %P[4]%x1' ^) ^
             ^( -clone 0 -process 'deintegim window %P[5]%x1' ^)
  set PROCS2=^( -clone 0 -process 'deintegim window 1x%P[4]%' ^) ^
             ^( -clone 0 -process 'deintegim window 1x%P[5]%' ^)
  set sPOLY=-poly "%W[1]%,1 %W[2]%,1 %W[3]%,1 %W[4]%,1 %W[5]%,1"
)

:: Next assumes k=4
:: Maybe we should add 8 transformations. This might explain the factor of 2 in the weights.

goto skip1

%IM7DEV%magick ^
  %INFILE% ^
  -process 'integim' ^
  ( -clone 0 -process 'deintegim window %P[1]%x1' -evaluate Multiply %W[1]% ) ^
  ( -clone 0 -process 'deintegim window %P[2]%x1' -evaluate Multiply %W[2]% ) ^
  ( -clone 0 -process 'deintegim window %P[3]%x1' -evaluate Multiply %W[3]% ) ^
  ( -clone 0 -process 'deintegim window %P[4]%x1' -evaluate Multiply %W[4]% ) ^
  -delete 0 ^
  -evaluate-sequence Add ^
+write i.png ^
  -process 'integim' ^
  ( -clone 0 -process 'deintegim window 1x%P[1]%' -evaluate Multiply %W[1]% ) ^
  ( -clone 0 -process 'deintegim window 1x%P[2]%' -evaluate Multiply %W[2]% ) ^
  ( -clone 0 -process 'deintegim window 1x%P[3]%' -evaluate Multiply %W[3]% ) ^
  ( -clone 0 -process 'deintegim window 1x%P[4]%' -evaluate Multiply %W[4]% ) ^
  -delete 0 ^
  +write d.png ^
  -evaluate-sequence Add ^
-auto-level ^
  %OUTFILE%

:skip1

:: Previous, but using "-poly".
%IM7DEV%magick ^
  %INFILE% ^
  -process 'integim' ^
  ( -clone 0 -process 'deintegim window %P[1]%x1' ) ^
  ( -clone 0 -process 'deintegim window %P[2]%x1' ) ^
  ( -clone 0 -process 'deintegim window %P[3]%x1' ) ^
  %PROCS1% ^
  -delete 0 ^
  %sPOLY% ^
  -process 'integim' ^
  ( -clone 0 -process 'deintegim window 1x%P[1]%' ) ^
  ( -clone 0 -process 'deintegim window 1x%P[2]%' ) ^
  ( -clone 0 -process 'deintegim window 1x%P[3]%' ) ^
  %PROCS2% ^
  -delete 0 ^
  %sPOLY% ^
-auto-level ^
  %OUTFILE%


goto skip2
:: Next looks blocky, horrible.
::
%IM7DEV%magick ^
  %INFILE% ^
  -process 'integim' ^
  ( -clone 0 -process 'deintegim window %P[1]%x1' -evaluate Multiply %Wadj[1]% ) ^
  ( -clone 0 -process 'deintegim window %P[2]%x1' -evaluate Multiply %Wadj[2]% ) ^
  ( -clone 0 -process 'deintegim window %P[3]%x1' -evaluate Multiply %Wadj[3]% ) ^
  ( -clone 0 -process 'deintegim window %P[4]%x1' -evaluate Multiply %Wadj[4]% ) ^
  ( -clone 0 -process 'deintegim window 1x%P[1]%' -evaluate Multiply %Wadj[1]% ) ^
  ( -clone 0 -process 'deintegim window 1x%P[2]%' -evaluate Multiply %Wadj[2]% ) ^
  ( -clone 0 -process 'deintegim window 1x%P[3]%' -evaluate Multiply %Wadj[3]% ) ^
  ( -clone 0 -process 'deintegim window 1x%P[4]%' -evaluate Multiply %Wadj[4]% ) ^
  -delete 0 ^
  +write d.png ^
  -evaluate-sequence Add ^
  %OUTFILE%
:skip2

endlocal

exit /B 0

::----------------------
:: Subroutines

:vplus1
set /A vp1=%1+1
set Wp1=!Wi[%vp1%]!

iiGaussSlice.bat

rem Find the best radius, hence diameter, of a block for a slice of a Gaussian kernel,
rem as measured by RMSE.

rem %1 sigma
rem %2 Vtop
rem %3 Vbottom

@if "%3"=="" findstr /B "rem @rem" %~f0 & exit /B 1

@setlocal enabledelayedexpansion

@call echoOffSave

echo %0: [%1] [%2] [%3]

set sigma=%1
if "%sigma%"=="." set sigma=
if "%sigma%"=="" set sigma=10

set Vtop=%2
set Vbottom=%3

set TMPDIR=\temp\
set HalfGauss=%TMPDIR%%Prefix%_iigkc2b.miff
set CurveSlice=%TMPDIR%%Prefix%_iicrvsl.miff


for /F "usebackq" %%L in (`%IMG7%magick ^
  -size %%[fx:%sigma%*10]x1 xc:Black ^
  -size 1x1 xc:White ^
  ^( -clone 0 ^) ^
  +append +repage ^
  -colorspace Gray ^
  -define "quantum:format=floating-point" -depth 32 ^
  -morphology Convolve Blur:0x%sigma% ^
  -alpha off ^
  -precision 15 ^
  -bordercolor Black -border 1 ^
  -trim +repage ^
  -format "WholeArea=%%[fx:mean*w]\n" +write info: ^
  -gravity East -crop %%[fx:^(w+1^)/2]x1+0+0 +repage ^
  +write %HalfGauss% ^
  -format "radP1=%%w\nMaxV=%%[fx:p[0,0]]\nHalfArea=%%[fx:mean*w]" ^
  info:`) do set %%L

if "%MaxV%"=="" (
  echo %0: MaxV not set
  exit /B 1
)

echo %0: radP1=%radP1% MaxV=%MaxV% Vtop=%Vtop% HalfArea=%HalfArea% WholeArea=%WholeArea%

set ReqWt=%Vtop%

for /F "usebackq" %%L in (`%IMG7%magick ^
  xc: ^
  -precision 15 ^
  -format "ThisWt=%%[fx:%Vtop%-!Vbottom!]\n" ^
  info:`) do set %%L


echo ThisWt=%ThisWt%

for /F "usebackq" %%L in (`%IMG7%magick ^
  xc: ^
  -precision 15 ^
  -format "ctop=%%[fx:!Vtop!>%MaxV%?!Vtop!:%MaxV%]" ^
  info:`) do set %%L

for /F "usebackq" %%L in (`%IMG7%magick ^
  %HalfGauss% ^
  -precision 15 ^
  -alpha off ^
  -fx "u>!Vtop!?%%[fx:!Vtop!-!Vbottom!]:u-!Vbottom!" ^
  -clamp ^
  -define "quantum:format=floating-point" -depth 32 ^
  +write %CurveSlice% ^
  -format "CurveSemiArea=%%[fx:mean*w]\n" +write info: ^
  -format "Xoffs=%%[fx:floor(mean*w/!ThisWt!-0.5)]\n" +write info: ^
  -format "BlockSemiArea=%%[fx:!ThisWt!*floor(mean*w/!ThisWt!-0.5)+!ThisWt!]\n" +write info: ^
  NULL:`) do set %%L


echo %0: ctop=!ctop! Vtop=!Vtop! Vbottom=!VBottom! CurveSemiArea=!CurveSemiArea! BlockSemiArea=!BlockSemiArea! ThisWt=!ThisWt! Xoffs=!Xoffs!


rem Some pixels from the left of the curve and block are equal,
rem so the boundary cannot occur there.

for /F "usebackq" %%L in (`%IMG7%magick ^
  %CurveSlice% ^
  -auto-level ^
  -fill Black +opaque White ^
  -format "StartRad=%%[fx:mean*w]\n" +write info: ^
  NULL:`) do set %%L

echo %0: StartRad=%StartRad% radP1=%radP1%

set BestComp=2.0
set Xoffs=0
set GoodEnuf=0

for /L %%I in (%StartRad%,1,%radP1%) do if !GoodEnuf!==0 (

  for /F "usebackq" %%L in (`%IMG7%magick ^
    -size %%Ix1 xc:White ^
    -evaluate Multiply %ThisWt% ^
    -background Black -extent %radP1% ^
    %CurveSlice% ^
    -metric RMSE ^
    -precision 15 ^
    -format %%[distortion] ^
    -compare ^
    info:`) do set comp=%%L

  for /F "usebackq" %%L in (`%IMG7%magick ^
    xc: ^
    -format "IsBetter=%%[fx:!BestComp!>!comp!?1:0]" ^
    info:`) do set %%L

  if !IsBetter!==1 (
    set BestComp=!comp!
    set Xoffs=%%I
  ) else (
    set GoodEnuf=1
  )
  rem echo %%I !comp! IsBetter=!IsBetter! GoodEnuf=!GoodEnuf!
)

set /A RatBest=%Xoffs%

set /A DatBest=%RatBest%*2+1

echo BestComp=%BestComp% RatBest=%RatBest% DatBest=%DatBest%

rem In next "+1" because we need to include the central pixel.

for /F "usebackq" %%L in (`%IMG7%magick ^
  xc: ^
  -precision 15 ^
  -alpha off ^
  -format "BlockSemiArea=%%[fx:%ThisWt%*(%RatBest%+1)]\n" ^
  info:`) do set %%L

for /F "usebackq" %%L in (`%IMG7%magick ^
  xc: ^
  -precision 15 ^
  -alpha off ^
  -format "AreaErr=%%[fx:abs(%CurveSemiArea%-%BlockSemiArea%)/%CurveSemiArea%]\n" ^
  info:`) do set %%L

echo CurveSemiArea=%CurveSemiArea% BlockSemiArea=%BlockSemiArea% AreaErr=%AreaErr% DatBest=%DatBest%

call echoRestore

@endlocal & set DatBest=%DatBest%

iiGaussKcalc.bat

rem What multiple different deintegrations of an integral image is closest to a 1D Gaussian blur?
rem Stacked Integral Images (SII) by the mid-edge method.
rem This is quite slow.

rem %1 Sigma [default 10]
rem %2 K, number of blocks [default 3]
rem %3 Prefix for output variables [default iigk_]

rem We can't "setlocal".
rem @setlocal enabledelayedexpansion

@call echoOffSave

set sigma=%1
if "%sigma%"=="." set sigma=
if "%sigma%"=="" set sigma=10

set K=%2
if "%K%"=="." set K=
if "%K%"=="" set K=3

set Prefix=%3
if "%Prefix%"=="." set Prefix=
if "%Prefix%"=="" set Prefix=iigk_

echo %0: Calculating kernel

set CNT=0

for /F "usebackq tokens=1,*" %%A in (`%IMG7%magick ^
  xc: ^
  -precision 15 ^
  -define "morphology:showKernel=1" ^
  -morphology Convolve:0 Blur:0x%sigma% ^
  info: 2^>^&1 ^|grep "0:" ^|head -n 1 ^| tr "  " \r\n `) do (
  if not "%%A"=="0:" (
    set KNL_V[!CNT!]=%%A
    set /A CNT+=1
  )
)

rem %CNT% should be an odd number

rem set KNL_V[


echo %0: CNT=%CNT%

set KNL_V[!CNT!]=0

set /A nLast=%CNT%-1
set /A nMid=%CNT%/2

echo %0: nMid=%nMid% nLast=%nLast%

set /A Km1=%K%-1

set Vtop=!KNL_V[%nMid%]!

set ReqWt=%Vtop%/%K%

if %K% gtr 1 set ReqWt=%ReqWt% + 0.25*%ReqWt%/%Km1%

for /F "usebackq" %%L in (`%IMG7%magick ^
  xc: ^
  -precision 15 ^
  -format "ReqWt=%%[fx:%ReqWt%]" ^
  info:`) do set %%L

echo %0: ReqWt=%ReqWt%

set Area=0

set startJ=%nMid%

for /L %%I in (0,1,%Km1%) do (
  if %%I==%Km1% (
    set ThisWt=!VTop!
    set Vbottom=0
  ) else (
    set ThisWt=%ReqWt%
    for /F "usebackq" %%L in (`%IMG7%magick ^
      xc: ^
      -precision 15 ^
      -format "ThisWt=%%[fx:%ReqWt%]\nVbottom=%%[fx:!Vtop!-^(%ReqWt%^)]\n" ^
      info:`) do set %%L
  )

  for /F "usebackq" %%L in (`%IMG7%magick ^
      xc: ^
      -precision 15 ^
      -format "Vmid=%%[fx:^(!Vtop!+!Vbottom!^)/2]" ^
      info:`) do set %%L

  @rem Find the value closest to Vmid. This is not fast.

  set MinDiff=999
  set JatMin=0
  set Done=0

  for /L %%J in (!startJ!,1,%nLast%) do (
    if !Done!==0 (
      for /F "usebackq" %%L in (`%IMG7%magick ^
        xc: ^
        -precision 15 ^
        -format "Diff=%%[fx:abs^(!Vmid!-!KNL_V[%%J]!^)]" ^
        info:`) do set %%L
      for /F "usebackq" %%L in (`%IMG7%magick ^
        xc: ^
        -precision 15 ^
        -format "Better=%%[fx:!Diff!<!MinDiff!?1:0]" ^
        info:`) do set %%L
      if !Better!==1 (
        set MinDiff=!Diff!
        set JatMin=%%J
      ) else (
        set Done=1
      )
      rem echo J=%%J  Diff=!Diff! Better=!Better! MinDiff=!MinDiff!
    )
  )

  echo %%I Vtop=!Vtop! Vbottom=!VBottom! Vmid=!Vmid! JatMin=!JatMin!

  set startJ=!JatMin!

  set Vtop=!Vbottom!

  set /A DatMin=^(!JatMin!-%nMid%^)*2+1

  for /F "usebackq" %%L in (`%IMG7%magick ^
      xc: ^
      -precision 15 ^
      -format "Area=%%[fx:!Area!+!DatMin!*!ThisWt!]" ^
      info:`) do set %%L

  set %Prefix%Diam[%%I]=!DatMin!
  set %Prefix%Weight[%%I]=!ThisWt!

  echo %0: DatMin=!DatMin! ThisWt=!ThisWt!
)

echo Area=%Area%

@rem Quantisation means the area is unlikely to be exactly 1.0.
@rem Adjust all the weights.

echo %0: Adjusting weights for area.

for /L %%I in (0,1,%Km1%) do (
  for /F "usebackq" %%L in (`%IMG7%magick ^
      xc: ^
      -precision 15 ^
      -format "%Prefix%Weight[%%I]=%%[fx:!%Prefix%Weight[%%I]!/%Area%]" ^
      info:`) do set %%L
)

set %Prefix%K=%K%
set %Prefix%Sigma=%Sigma%

call echoRestore

@endlocal

iiGaussKcalc_4.bat

rem Like iiGaussCalc.bat, but finds 50% of area, processing by images.

rem What multiple different deintegrations of an integral image is closest to a 1D Gaussian blur?

rem %1 Sigma [default 10]
rem %2 K, number of blocks [default 3]
rem %3 Prefix for output variables [default iigk_]
rem %4 Limit area error (0 to 1) [default 0.0001]

@rem This works from the lowest (widest) slice upwards.


@setlocal enabledelayedexpansion

@call echoOffSave


set sigma=%1
if "%sigma%"=="." set sigma=
if "%sigma%"=="" set sigma=10

set K=%2
if "%K%"=="." set K=
if "%K%"=="" set K=3

set Prefix=%3
if "%Prefix%"=="." set Prefix=
if "%Prefix%"=="" set Prefix=iigk_

set LimAreaErr=%4
if "%LimAreaErr%"=="." set LimAreaErr=
if "%LimAreaErr%"=="" set LimAreaErr=0.0001

set TMPDIR=\temp\
set HalfGauss=%TMPDIR%%Prefix%_iigkc2b.miff

rem echo %0: Calculating kernel

@rem Vtop etc are on a scale of 0.0 to 1.0.

set /A Km1=%K%-1

set Vtop=

for /F "usebackq" %%L in (`%IMG7%magick ^
  -size %%[fx:%sigma%*10]x1 xc:Black ^
  -size 1x1 xc:White ^
  ^( -clone 0 ^) ^
  +append +repage ^
  -colorspace Gray ^
  -define "quantum:format=floating-point" -depth 32 ^
  -morphology Convolve Blur:0x%sigma% ^
  -alpha off ^
  -precision 15 ^
  -bordercolor Black -border 1 ^
  -trim +repage ^
  -format "WholeArea=%%[fx:mean*w]\n" +write info: ^
  -gravity East -crop %%[fx:^(w+1^)/2]x1+0+0 +repage ^
  +write %HalfGauss% ^
  -format "radP1=%%w\nMaxV=%%[fx:p[0,0]]\nHalfArea=%%[fx:mean*w]" ^
  info:`) do set %%L

if "%MaxV%"=="" (
  echo %0: MaxV not set
  exit /B 1
)
set Vtop=%MaxV%

echo %0: radP1=%radP1% Vtop=%Vtop% HalfArea=%HalfArea% WholeArea=%WholeArea%

set ReqWt=%Vtop%/%K%

if %K%==1 (
  set ReqWtBot=%Vtop%
  set ReqWtOth=%Vtop%
) else (
  set ReqWtBot=%Vtop%/%K%/1.5
  set "ReqWtOth=(%Vtop%-!ReqWtBot!)/%Km1%"
)

rem if %K% gtr 1 set ReqWt=%ReqWt% + 0.25*%ReqWt%/%Km1%

for /F "usebackq" %%L in (`%IMG7%magick ^
  xc: ^
  -precision 15 ^
  -format "ReqWt=%%[fx:%ReqWt%]\n" ^
  info:`) do set %%L

for /F "usebackq" %%L in (`%IMG7%magick ^
  xc: ^
  -precision 15 ^
  -format "ReqWtOth=%%[fx:%ReqWtOth%]\n" ^
  info:`) do set %%L

for /F "usebackq" %%L in (`%IMG7%magick ^
  xc: ^
  -precision 15 ^
  -format "ReqWtBot=%%[fx:%ReqWtBot%]\n" ^
  info:`) do set %%L

echo %0: Vtop=%Vtop% ReqWt=%ReqWt% ReqWtOth=%ReqWtOth% ReqWtBot=%ReqWtBot%

set Area=0

for /L %%I in (%Km1%,-1,0) do (
  if %%I==%Km1% (
    set Vtop=%ReqWtBot%
    set ThisWt=!VTop!
    set Vbottom=0
  ) else if %%I==0 (
    for /F "usebackq" %%L in (`%IMG7%magick ^
      xc: ^
      -precision 15 ^
      -format "Vtop=%%[fx:%MaxV%]\n" -write info: ^
      -format "ThisWt=%%[fx:%MaxV%-!Vbottom!]\n" ^
      info:`) do set %%L
  ) else (
    set ThisWt=%ReqWtOth%
    for /F "usebackq" %%L in (`%IMG7%magick ^
      xc: ^
      -precision 15 ^
      -format "Vtop=%%[fx:!Vbottom!+%ReqWtOth%]\n" ^
      info:`) do set %%L
  )

  echo %%I Vtop=!Vtop! Vbottom=!Vbottom! ThisWt=!ThisWt!

  @rem Ignoring values above Vtop and below Vbottom, just for the RHS, where is the area half of total?.

  set Xoffs=

  set GoodEnuf=0
  set nIter=0

  for /L %%J in (0,1,99) do if !GoodEnuf!==0 (

    if %%I==0 (
      for /F "usebackq" %%L in (`%IMG7%magick ^
        xc: ^
        -precision 15 ^
        -format "ctop=%%[fx:!Vtop!>%MaxV%?!Vtop!:%MaxV%]" ^
        info:`) do set %%L

    ) else (
      set ctop=!Vtop!
    )

    for /F "usebackq" %%L in (`%IMG7%magick ^
      %HalfGauss% ^
      -precision 15 ^
      -alpha off ^
      -fx "u>!ctop!?%%[fx:!ctop!-!Vbottom!]:u-!Vbottom!" ^
      -clamp ^
      -format "CurveSemiArea=%%[fx:mean*w]\n" +write info: ^
      -format "Xoffs=%%[fx:floor(mean*w/!ThisWt!-0.5)]\n" +write info: ^
      -format "BlockSemiArea=%%[fx:!ThisWt!*floor(mean*w/!ThisWt!-0.5)+!ThisWt!]\n" +write info: ^
      NULL:`) do set %%L

rem       -evaluate Subtract %%[fx:!Vbottom!*QuantumRange] -clamp

    if "!Xoffs!"=="" (
      echo %0: Xoffs not set. I=%%I
      exit /B 1
    )
    rem echo %0: %%I Vtop=!Vtop! Vbottom=!VBottom! CurveSemiArea=!CurveSemiArea! BlockSemiArea=!BlockSemiArea! ThisWt=!ThisWt! Xoffs=!Xoffs!

    rem The Xoffs must be an integer, so the area will have a quantisation error.
    rem Adjust this weight to compensate.
    rem We start from bottom slice.

    for /F "usebackq" %%L in (`%IMG7%magick ^
      xc: ^
      -format "GoodEnuf=%%[fx:abs(!CurveSemiArea!-!BlockSemiArea!)/!CurveSemiArea!>%LimAreaErr%?0:1]\n" ^
      info:`) do set %%L

    set /A nIter+=1
    set DontTrim=0

    if !GoodEnuf!==0 (

      for /F "usebackq" %%L in (`%IMG7%magick ^
        xc: ^
        -precision 15 ^
        -format "ThisWt=%%[fx:!ThisWt!*!CurveSemiArea!/!BlockSemiArea!]\n" ^
        info:`) do set %%L

      for /F "usebackq" %%L in (`%IMG7%magick ^
        xc: ^
        -precision 15 ^
        -format "Vtop=%%[fx:!Vbottom!+!ThisWt!]\n" +write info: ^
        info:`) do set %%L

rem      if %%I==0 (
rem        for /F "usebackq" %%L in (`%IMG7%magick ^
rem          xc: ^
rem          -precision 15 ^
rem          -format "DontTrim=%%[fx:!tmpVtop!<%MaxV%?1:0]\n" +write info: ^
rem          info:`) do set %%L
rem      )
rem
rem      if !DontTrim!==1 (
rem        set GoodEnuf=1
rem      ) else (
rem        set ThisWt=!tmpThisWt!
rem        set Vtop=!tmpVtop!
rem      )
rem
rem      set ThisWt=!tmpThisWt!
rem      set Vtop=!tmpVtop!

      rem echo %0: new ThisWt=!ThisWt! Vtop=!Vtop! GoodEnuf=!GoodEnuf!

      rem The new weight and Vtop means the slice covers a different part of the curve,
      rem so go round again.
    )
  )
  rem End J loop

  set Vbottom=!Vtop!

  set /A DatMin=Xoffs*2+1

  for /F "usebackq" %%L in (`%IMG7%magick ^
      xc: ^
      -precision 15 ^
      -format "Area%%I=%%[fx:!DatMin!*!ThisWt!]" ^
      info:`) do set %%L

  echo %0: DatMin=!DatMin! ThisWt=!ThisWt! CurveSemiArea=!CurveSemiArea! BlockSemiArea=!BlockSemiArea! Area%%I=!Area%%I! GoodEnuf=!GoodEnuf! nIter=!nIter!

  for /F "usebackq" %%L in (`%IMG7%magick ^
      xc: ^
      -precision 15 ^
      -format "Area=%%[fx:!Area!+!DatMin!*!ThisWt!]" ^
      info:`) do set %%L

  set %Prefix%Diam[%%I]=!DatMin!
  set %Prefix%Weight[%%I]=!ThisWt!

  rem echo %0: DatMin=!DatMin! ThisWt=!ThisWt!
)

if ERRORLEVEL 1 exit /B 1

echo %0 Area=%Area%

set %Prefix%K=%K%
set %Prefix%Sigma=%sigma%

@rem Quantisation means the total area is unlikely to be exactly 1.0.
@rem Adjust all the weights.

goto skipAdj

set %Prefix%

echo %0: Adjusting weights for area.

for /L %%I in (0,1,%Km1%) do (
  for /F "usebackq" %%L in (`%IMG7%magick ^
      xc: ^
      -precision 15 ^
      -format "%Prefix%Weight[%%I]=%%[fx:!%Prefix%Weight[%%I]!/%Area%]" ^
      info:`) do set %%L
)

set %Prefix%

:skipAdj

set %Prefix% > %Prefix%_vars.lis

call echoRestore

@endlocal & set Prefix=%Prefix%

@for /F %%L in (%Prefix%_vars.lis) do @set %%L

iiGaussK.bat

rem Blurs an image using the Stacked Integral Images (SII) method
rem from parameters calculated by iiGaussKcalc.bat.

rem %1 input image
rem %2 output image
rem %3 Prefix for variables already calculated [default iigk_]

@setlocal enabledelayedexpansion

@call echoOffSave

set INFILE=%1
set OUTFILE=%2

set Prefix=%3
if "%Prefix%"=="." set Prefix=
if "%Prefix%"=="" set Prefix=iigk_


echo %0: Prefix=%Prefix%

set %Prefix%

if "!%Prefix%K!" == "" (
  if not exist %Prefix%_vars.lis (
    echo Can't find file %Prefix%_vars.lis
    exit /B 1
  )
  echo %0: Reading env vars from %Prefix%_vars.lis
  @for /F %%L in (%Prefix%_vars.lis) do @set %%L

)

if "!%Prefix%K!" == "" (
  echo %0: %Prefix%K is blank.
  exit /B 1
)

set TMPDIR=\temp

set fDeinteg1=%TMPDIR%\%Prefix%1.txt
set fDeinteg2=%TMPDIR%\%Prefix%2.txt
set fPoly=%TMPDIR%\%Prefix%poly.txt
set fScr=%TMPDIR%\%Prefix%.scr

if "!%Prefix%K!" == "" (
  echo %0: iiGaussKcalc needs to be successfully called.
  exit /B 1
)

set K=!%Prefix%K!
set /A Km1=%K%-1

(
  for /L %%I in (0,1,%Km1%) do echo ^( -clone 0 -process 'deintegim window !%Prefix%Diam[%%I]!x1 adjedges' ^)
) >%fDeinteg1%

(
  for /L %%I in (0,1,%Km1%) do echo ^( -clone 0 -process 'deintegim window 1x!%Prefix%Diam[%%I]! adjedges' ^)
) >%fDeinteg2%

(
  for /L %%I in (0,1,%Km1%) do echo !%Prefix%Weight[%%I]!,1,
) >%fPoly%


echo %0: Writing %fScr%
(
  echo -process integim
  type %fDeinteg1%
  echo -delete 0
  echo -poly @"%fPoly%"
  echo -process integim
  type %fDeinteg2%
  echo -delete 0
  echo -poly @"%fPoly%"
  echo -write %OUTFILE%
  echo -exit
) >%fScr%

echo %0: Executing IM: INFILE=%INFILE%

type %fScr%

%IM7DEVFX3%magick %INFILE% -script %fScr%

call echoRestore

@endlocal

iiGaussKgraph.bat

rem From data calculated by iiGaussKcalc,
rem make a graph.

rem %1 Prefix for variables already calculated [default iigk_]
rem %2 output filename

@setlocal enabledelayedexpansion

@call echoOffSave

set Prefix=%1
if "%Prefix%"=="." set Prefix=
if "%Prefix%"=="" set Prefix=iigk_

set OUTFILE=%2
if "%OUTFILE%"=="." set OUTFILE=
if "%OUTFILE%"=="" set OUTFILE=iigkg.png

if "!%Prefix%K!" == "" (
  if not exist %Prefix%_vars.lis (
    echo Can't find file %Prefix%_vars.lis
    exit /B 1
  )
  echo %0: Reading env vars from %Prefix%_vars.lis
  @for /F %%L in (%Prefix%_vars.lis) do @set %%L

)

if "!%Prefix%K!" == "" (
  echo %0: %Prefix%K is blank.
  exit /B 1
)

set %Prefix%

set TMPDIR=\temp\
set TMPFILE=%TMPDIR%%Prefix%iigkg.tiff
set TMPOUT=%TMPDIR%%Prefix%iigkg_out.tiff

set K=!%Prefix%K!
set /A Km1=%K%-1

set sig=!%Prefix%Sigma!

for /F "usebackq" %%L in (`%IMG7%magick ^
  -size %%[fx:%sig%*10]x1 xc:Black ^
  -size 1x1 xc:White ^
  ^( -clone 0 ^) ^
  +append +repage ^
  -morphology Convolve Blur:0x%sig% ^
  -bordercolor Black -border 1 ^
  -trim +repage ^
  -set option:Xscale "%%[fx:w<100?6:w<200?3:1]" ^
  -scale "%%[fx:Xscale*100]x100%%" ^
  -precision 15 ^
  -set option:Vmax %%[fx:maxima*1.1] ^
  -evaluate Divide %%[fx:Vmax] ^
  -format "Xscale=%%[Xscale]\nVmax=%%[Vmax]\n" +write info: ^
  -define "quantum:format=floating-point" ^
  -depth 32 ^
  %TMPFILE%`) do set %%L

if "%Vmax%"=="" (
  echo %0: Vmax is not set. sig=%sig%
  exit /B 1
)

rem But Xscale is messed up by graphLineCol

echo %0: Vmax=%Vmax%

call %PICTBAT%graphLineCol %TMPFILE% . . 0 %TMPOUT%

for /F "usebackq" %%L in (`%IMG7%magick ^
  %TMPOUT% ^
  -gravity NorthWest -fill White -pointsize 15 -annotate +10+10 "sigma=%sig%\nK=%K%" ^
  -flip ^
  +write %TMPOUT% ^
  -format "Xscale=%%[fx:%Xscale%*%glcXscale%]\nYscale=%%[fx:h/%Vmax%]\n" ^
  info:`) do set %%L

echo %0: Xscale=%Xscale% Yscale=%Yscale%

set Cols[0]=#08f8
set Cols[1]=#4f48
set Cols[2]=#f448

set Vbottom=0

for /L %%I in (%Km1%,-1,0) do (
  set ThisWt=!%Prefix%Weight[%%I]!
  set Vtop=!Vbottom!+!ThisWt!
  set Diam=!%Prefix%Diam[%%I]!
  set /A Rad=^(!Diam!-1^)/2

  set /A ColNdx=%%I%%3
  set Col=Cols[!ColNdx!]
  set ColV=!Col!

  rem FIXME: we might adjust rectangle top or bottom to avoid overlaps.

  call :getCol %%I

  echo %0: I=%%I ThisWt=!ThisWt! Diam=!Diam! Rad=!Rad! Vbottom=!Vbottom! ColNdx=!ColNdx! Col=!Col! ColV=!ColV!

  for /F "usebackq" %%L in (`%IMG7%magick ^
    %TMPOUT% ^
    -precision 15 ^
    -set option:Vtop %%[fx:!Vtop!] ^
    -stroke None -fill !Col! ^
    -draw "rectangle %%[fx:w/2-%Xscale%*!Rad!],%%[fx:!Vbottom!*%Yscale%],%%[fx:w/2+%Xscale%*!Rad!],%%[fx:Vtop*%Yscale%]" ^
    -format "Vtop=%%[Vtop]" +write info: ^
    %TMPOUT%`) do set %%L

  set Vbottom=!Vtop!
)

%IMG7%magick ^
  %TMPOUT% ^
  -flip ^
  -bordercolor gray(50%%) -border 5 ^
  %OUTFILE%

call echoRestore

@endlocal & set iigkgOUTFILE=%OUTFILE%

exit /B 0

::-------------------------------------------------------------
:: Subroutines

:getCol
  set /A ColNdx=%1%%3
  set Col=!Cols[%ColNdx%]!

  exit /B 0

iiGaussChkRmse.bat

rem %1 Input prefix for variables already calculated [default iigk_]
rem %2 Output prefix for new values [default no output]

@setlocal enabledelayedexpansion

@call echoOffSave

set Prefix=%1
if "%Prefix%"=="." set Prefix=
if "%Prefix%"=="" set Prefix=iigk_

set OutPrefix=%2
if "%OutPrefix%"=="." set OutPrefix=


echo %0: Prefix=%Prefix%

set %Prefix%

if "!%Prefix%K!" == "" (
  if not exist %Prefix%_vars.lis (
    echo Can't find file %Prefix%_vars.lis
    exit /B 1
  )
  echo %0: Reading env vars from %Prefix%_vars.lis
  @for /F %%L in (%Prefix%_vars.lis) do @set %%L

)

if "!%Prefix%K!" == "" (
  echo %0: %Prefix%K is blank.
  exit /B 1
)

set sigma=!%Prefix%Sigma!
set K=!%Prefix%K!
set /A Km1=%K%-1

set Vbottom=0

if not "%OutPrefix%"=="" (
  set %OutPrefix%K=%K%
  set %OutPrefix%Sigma=%sigma%
)

for /L %%I in (%Km1%,-1,0) do (
  rem call :WtNext %%I
  rem set /A Ip1=%%I+1
  rem set Vtop=!%Prefix%Weight[%%I]!
  rem set VBottom=!Vnext!

  set ThisWt=!%Prefix%Weight[%%I]!
  set %OutPrefix%Weight[%%I]=!ThisWt!

  for /F "usebackq" %%L in (`%IMG7%magick ^
    xc: ^
    -format "Vtop=%%[fx:!Vbottom!+!ThisWt!]\n" ^
    info:`) do set %%L


  echo %%I sigma=%sigma% Vtop=!Vtop! Vbottom=!Vbottom!

  call %PICTBAT%iiGaussSlice %sigma% !Vtop! !Vbottom!

  set %OutPrefix%Diam[%%I]=!DatBest!

  set Vbottom=!Vtop!
)

if not "%OutPrefix%"=="" (
  set %OutPrefix%

  set %OutPrefix% > %OutPrefix%_vars.lis

)


call echoRestore

@endlocal & set OutPrefix=%OutPrefix%

@for /F %%L in (%OutPrefix%_vars.lis) do @set %%L

@exit /B 0

::----------------------------------------------
:: Subroutine

:WtNext
set /A Ip1=%1+1
set WtNext=!%Prefix%Weight[%Ip1%]!
exit /B 0

iiGaussSlice.bat

rem Find the best radius, hence diameter, of a block for a slice of a Gaussian kernel,
rem as measured by RMSE.

rem %1 sigma
rem %2 Vtop
rem %3 Vbottom

@if "%3"=="" findstr /B "rem @rem" %~f0 & exit /B 1

@setlocal enabledelayedexpansion

@call echoOffSave

echo %0: [%1] [%2] [%3]

set sigma=%1
if "%sigma%"=="." set sigma=
if "%sigma%"=="" set sigma=10

set Vtop=%2
set Vbottom=%3

set TMPDIR=\temp\
set HalfGauss=%TMPDIR%%Prefix%_iigkc2b.miff
set CurveSlice=%TMPDIR%%Prefix%_iicrvsl.miff


for /F "usebackq" %%L in (`%IMG7%magick ^
  -size %%[fx:%sigma%*10]x1 xc:Black ^
  -size 1x1 xc:White ^
  ^( -clone 0 ^) ^
  +append +repage ^
  -colorspace Gray ^
  -define "quantum:format=floating-point" -depth 32 ^
  -morphology Convolve Blur:0x%sigma% ^
  -alpha off ^
  -precision 15 ^
  -bordercolor Black -border 1 ^
  -trim +repage ^
  -format "WholeArea=%%[fx:mean*w]\n" +write info: ^
  -gravity East -crop %%[fx:^(w+1^)/2]x1+0+0 +repage ^
  +write %HalfGauss% ^
  -format "radP1=%%w\nMaxV=%%[fx:p[0,0]]\nHalfArea=%%[fx:mean*w]" ^
  info:`) do set %%L

if "%MaxV%"=="" (
  echo %0: MaxV not set
  exit /B 1
)

echo %0: radP1=%radP1% MaxV=%MaxV% Vtop=%Vtop% HalfArea=%HalfArea% WholeArea=%WholeArea%

set ReqWt=%Vtop%

for /F "usebackq" %%L in (`%IMG7%magick ^
  xc: ^
  -precision 15 ^
  -format "ThisWt=%%[fx:%Vtop%-!Vbottom!]\n" ^
  info:`) do set %%L


echo ThisWt=%ThisWt%

for /F "usebackq" %%L in (`%IMG7%magick ^
  xc: ^
  -precision 15 ^
  -format "ctop=%%[fx:!Vtop!>%MaxV%?!Vtop!:%MaxV%]" ^
  info:`) do set %%L

for /F "usebackq" %%L in (`%IMG7%magick ^
  %HalfGauss% ^
  -precision 15 ^
  -alpha off ^
  -fx "u>!Vtop!?%%[fx:!Vtop!-!Vbottom!]:u-!Vbottom!" ^
  -clamp ^
  -define "quantum:format=floating-point" -depth 32 ^
  +write %CurveSlice% ^
  -format "CurveSemiArea=%%[fx:mean*w]\n" +write info: ^
  -format "Xoffs=%%[fx:floor(mean*w/!ThisWt!-0.5)]\n" +write info: ^
  -format "BlockSemiArea=%%[fx:!ThisWt!*floor(mean*w/!ThisWt!-0.5)+!ThisWt!]\n" +write info: ^
  NULL:`) do set %%L


echo %0: ctop=!ctop! Vtop=!Vtop! Vbottom=!VBottom! CurveSemiArea=!CurveSemiArea! BlockSemiArea=!BlockSemiArea! ThisWt=!ThisWt! Xoffs=!Xoffs!


rem Some pixels from the left of the curve and block are equal,
rem so the boundary cannot occur there.

for /F "usebackq" %%L in (`%IMG7%magick ^
  %CurveSlice% ^
  -auto-level ^
  -fill Black +opaque White ^
  -format "StartRad=%%[fx:mean*w]\n" +write info: ^
  NULL:`) do set %%L

echo %0: StartRad=%StartRad% radP1=%radP1%

set BestComp=2.0
set Xoffs=0
set GoodEnuf=0

for /L %%I in (%StartRad%,1,%radP1%) do if !GoodEnuf!==0 (

  for /F "usebackq" %%L in (`%IMG7%magick ^
    -size %%Ix1 xc:White ^
    -evaluate Multiply %ThisWt% ^
    -background Black -extent %radP1% ^
    %CurveSlice% ^
    -metric RMSE ^
    -precision 15 ^
    -format %%[distortion] ^
    -compare ^
    info:`) do set comp=%%L

  for /F "usebackq" %%L in (`%IMG7%magick ^
    xc: ^
    -format "IsBetter=%%[fx:!BestComp!>!comp!?1:0]" ^
    info:`) do set %%L

  if !IsBetter!==1 (
    set BestComp=!comp!
    set Xoffs=%%I
  ) else (
    set GoodEnuf=1
  )
  rem echo %%I !comp! IsBetter=!IsBetter! GoodEnuf=!GoodEnuf!
)

set /A RatBest=%Xoffs%

set /A DatBest=%RatBest%*2+1

echo BestComp=%BestComp% RatBest=%RatBest% DatBest=%DatBest%

rem In next "+1" because we need to include the central pixel.

for /F "usebackq" %%L in (`%IMG7%magick ^
  xc: ^
  -precision 15 ^
  -alpha off ^
  -format "BlockSemiArea=%%[fx:%ThisWt%*(%RatBest%+1)]\n" ^
  info:`) do set %%L

for /F "usebackq" %%L in (`%IMG7%magick ^
  xc: ^
  -precision 15 ^
  -alpha off ^
  -format "AreaErr=%%[fx:abs(%CurveSemiArea%-%BlockSemiArea%)/%CurveSemiArea%]\n" ^
  info:`) do set %%L

echo CurveSemiArea=%CurveSemiArea% BlockSemiArea=%BlockSemiArea% AreaErr=%AreaErr% DatBest=%DatBest%

call echoRestore

@endlocal & set DatBest=%DatBest%

iiGaussOpn.bat

rem %1 Name of environment variable
rem %2 method: one of blur, morph, single2d, single1d, repeatedii, SIImid
rem %3 sigma
rem %4 K (number of slices)

@if "%1"=="" findstr /B "rem @rem" %~f0 & exit /B 1

@setlocal enabledelayedexpansion

rem @call echoOffSave

call %PICTBAT%setInOut %1 iigo

set ENVVAR=%1
if "%ENVVAR%"=="." set ENVVAR=
if "%ENVVAR%"=="" set ENVVAR=iigo_

set METH=%2
if "%METH%"=="." set METH=
if "%METH%"=="" set METH=blur

set SIGMA=%3
if "%SIGMA%"=="." set SIGMA=
if "%SIGMA%"=="" set SIGMA=10

set K=%4
if "%K%"=="." set K=
if "%K%"=="" set K=3

set sOPN=

if /I "%METH%"=="blur" (
  set sOPN=-blur 0x%SIGMA%

) else if /I "%METH%"=="morph" (
  set sOPN=-morphology Convolve "Blur:0x%SIGMA%;Blur:0x%SIGMA%,90"

) else if /I "%METH%"=="single2d" (
  for /F "usebackq" %%L in (`%IMG7%magick ^
    xc: ^
    -format "calcW=%%[fx:calcR=max(1,floor(%SIGMA%*1.63)); calcR-=(%SIGMA%>100)?floor(0.5+575*pow((%SIGMA%-100)/1900,1.409)):0; 2*calcR+1]" ^
    info:`) do set %%L
  set sOPN=-process integim -process 'deintegim window !calcW!x!calcW!'

) else if /I "%METH%"=="single1d" (
  for /F "usebackq" %%L in (`%IMG7%magick ^
    xc: ^
    -format "bestW=%%[fx:2 * floor (1.53 * %SIGMA%) + 1]" ^
    info:`) do set %%L
  set sOPN=-process integim -process 'deintegim window !calcW!x!calcW!'

) else if /I "%METH%"=="repeatedii" (
  for /F "usebackq" %%L in (`%IMG7%magick ^
    xc: ^
    -format "bestW=%%[fx:2 * floor (0.95 * %SIGMA%) + 1]" ^
    info:`) do set %%L
  set sOPN1=-process integim -process 'deintegim window !bestW!x1'
  set sOPN2=-process integim -process 'deintegim window 1x!bestW!'
  set sOPN=!sOPN1! !sOPN1! !sOPN1! !sOPN2! !sOPN2! !sOPN2! 

) else if /I "%METH%"=="SIImid" (
  call %PICTBAT%iiGaussKcalc %SIGMA% 3 iigk_

  set sOPN=-blur 0x%SIGMA%
)

echo %0: %sOPN%

call echoRestore

endlocal & set iigoOUTFILE=%OUTFILE%

%IMG7%magick -version 

Version: ImageMagick 7.1.0-49 Q16-HDRI x64 7a3f3f1:20220924 https://imagemagick.org
Copyright: (C) 1999 ImageMagick Studio LLC
License: https://imagemagick.org/script/license.php
Features: Cipher DPC HDRI OpenCL 
Delegates (built-in): bzlib cairo freetype gslib heic jng jp2 jpeg jxl lcms lqr lzma openexr pangocairo png ps raqm raw rsvg tiff webp xml zip zlib
Compiler: Visual Studio 2022 (193331630)

%IM7DEV%magick -version 

Version: ImageMagick 7.1.0-33 beta Q32-HDRI x86_64 64b5fe68a:20220501 https://imagemagick.org
Copyright: (C) 1999 ImageMagick Studio LLC
License: https://imagemagick.org/script/license.php
Features: Cipher DPC HDRI Modules OpenMP(4.5) 
Delegates (built-in): bzlib fpx jbig jpeg lcms ltdl lzma tiff webp x xml zip zlib
Compiler: gcc (11.3)

Page version v2.0 4-November-2022.

Page created 17-Mar-2023 03:21:37.

Copyright © 2023 Alan Gibson.