Sigma*

[emery]

The Imagick API doesn't have a default value for the standard deviation of any functions that take it, i.e. blur, and sharpen. The command line however does not require a stigma.

What is the default value so I don't have to require it in my wrapper?

Thank you

[anthony]

It is Sigma not Stigma -- I don't want bloody hands!

This problem has been noted many times before.

Basically the argument order of -blur {radius}x{sigma} makes 'sigma' optional when it isn't.
Also 'sigma' is the main control, not 'radius', and as such should be set 'first'. However it has always been like this and it is FAR to late to change it now.

However their is no GOOD default value you can use for sigma either!

[fmw42]

Actually if you want a linear blur rather than gaussian, you can use

radiusxsigma = radiusx65000

or any number for sigma that is about 10 times or more that of radius.


For gaussian blur,

radiusxsigma => 0xsigma will have a radius of about 3*sigma

[emery]

Yikes, forgive my misspelling =)

What would be the closest to photoshop? Also, photoshop has a "spread" argument, which I am emulating by blurring a stroke. Is there a way to get that directly from blur with radius and sigma?

[fmw42]

You need to explain more about what you are doing in Photoshop. Not everyone is familiar with its features. What and where do you find "spread" in Photoshop and what does it do?

[emery]

My mistake, I neglected to mention my intentions. I am actually going for an "outer glow" effect similar to that in photoshop.



The "size" option indicates the range of the blur and "spread" indicates how much of the size is actually blurred.

[fmw42]

Unfortunately I only have PS CS which does not seem to have that function. Take this image and process it and post a screen snap of the controls and the output from your processing. Perhaps a profile will tell me more what they are doing.

[emery]

Your version should have it. Under the "Layers" window pane there is an "fx" button which opens that Layer Styles window.




There is also a "Contour" setting which allows you to select separate waveform, and adjust the range:

[fmw42]

This is much more complex than anything in IM and certainly not just a -blur with radiusxsigma as you have a contour in PS in the controls that is not part of -blur. One may be able to reproduce something that matches the linear coutour with -blur or -gaussian-blur. But I would need to have access to the resulting images separately from your displayed images so that I don't have to crop them out and add further resampling. So if you could display the results separately from the control panel that would help. Getting the other contours may be possible with other IM processing but lets start with something simple. Also I have no idea what the technique "Softer" means nor what other options there are. Likewise the blend mode of "Screen". I don't know what the blending means nor what other choices you have.

You are trying to reproduce some complex operation involving color, blend, (opacity), (noise), technique, spread, size, contour, range, jitter -- that is 10 controls -- with something in IM that has only two controls, radius and sigma which may even be related, so perhaps as few as one control. In order to achieve any success, one needs to vary one item at a time.

So it would be best to keep contour at linear, pick one blend mode, keep opacity=100, noise=0, color=white, spread=0, range=50% and jitter=0. Then vary the size and lets see what you get.

Later we can try fixing all the others as above, keep the size fixed and then vary the spread to see what that does.

[emery]

It is worth noting that I am always working with vector graphics that are composited into a scalar image. So I have the ability to extend the perimeter of the shape with a stroke. I am able to reproduce something similar by blurring the stroke extended shape. It doesn't seem like a good technique though.

Let me also explain, that I'm merely trying to make my controls simpler for designers, since that is who will be using the end product, and most designers use photoshop. It's not important to have all of the same functionality, but I want to do everything I can to soften the learning curve for designers.

• "blend" is just the composite operator of the result, which I already handle.
• "noise" is just an after effect, already handled separately. They only have it there because the independent "noise" effect doesn't apply to the outer glow. You would have to collapse the layer first. Not a problem for me.
• "technique" is a minor insignificance in the algorithm, the result looks mostly the same.
• "jitter" is actually some other kind of noise, and also insignificant.
• "opacity" is needless because the color of the image to be blurred is already RGBA.

For now all I really care about is spread and size. The other controls are rarely ever used.


Here are those images:
size=10:


size=20:


size=40:


It does look like blurring doesn't it? If it was a special algorithm it would respect those corners.

[fmw42]

It looks like a circular blur (probably gaussian profile) from the center, then that blended with the square.
It could be that size is like radius and spread is like sigma on the Gaussian profile of the circularly symmetric blur.

Here is an example of generating a gaussian blur with circular shape and adding it onto the test square. The gaussian is generated with -blur from a white point in a black background and thus needs to be stretched to full dynamic range (as the more it is blurred the darker it will be). I used -linear-stretch 1x1 to achieve a stretch to full dynamic range. Then I multiply by 1.3 to enhance its "whiteness" and add to the test square.

test square:


Gaussian blur added to square:
convert testsq.png \
\( -size 128x128 xc:black -fill white \
-draw "point 64,64" -blur 0x30 -linear-stretch 1x1 -evaluate multiply 1.3 \) \
-compose plus -composite testsq_blur0x30_mult1p3.png



You can play with the radiusxsigma in blur. I used radius=0 so the radius will be generated from the sigma as approx. 3x30=90 (so should be similar to 90x30) and also the multiplier where I used 1.3

Also if radius and sigma don't give enough control, you can also apply a sigmoidal-contrast to the gaussian to make it more contrasty and change the profile. see -sigmoidal-contrast

other functions such as applying a sinusoid (see -evalute sine or -function sinusoid) could also apply the ripple effects that they use with other profiles.

[anthony]

I think fred's idea is the right one.

• Generate the blur by extracting the alpha channel and using a linear blur.
• Separately generate a contour as a gradient image (color and transparency).
• The first image is then used to do a -clut (color table lookup) of the contour image.
• overlay the original image.

That appears to be how it is all working.

However there is an alternative way to using a 'blur image'. That is by creating a 'distance from object type function. With that type of function the 'corners' of the shape would not 'stick out' from the glow on large blur radii. Unfortunately this type of function is not each to create, though it is in the 'morphology' methods of image handling, it is not available yet in IM, though it is a long term future proposal, as part of morphological operators.

[thevinn]

However there is an alternative way to using a 'blur image'. That is by creating a 'distance from object type function. With that type of function the 'corners' of the shape would not 'stick out' from the glow on large blur radii. Unfortunately this type of function is not each to create, though it is in the 'morphology' methods of image handling, it is not available yet in IM, though it is a long term future proposal, as part of morphological operators.

I believe this is the correct approach. It can be easily proven that Photohsop does not use a simple 8-bit mask to colourize the destination layer using the gradient colours: just make a blend with many stops, apply outer glow at full size (250px) with the blend, and sample the colour values. You will see more than 256 of them.

On the other hand, with a distance transform you can have a colour lookup table with any number of entries you want. I've managed to implement THREE separate, fast distance transform algorithms from various recent papers, here they are:

Code: Select all

/** Distance transform calculations.

    @ingroup vf_gui
*/
struct DistanceTransform
{
  //----------------------------------------------------------------------------

  /** Mask inclusion test functor.
  */
  struct WhiteTest
  {
    explicit WhiteTest (Pixels src) : m_src (src)
    {
    }

    inline bool operator () (int const x, int const y) const noexcept
    {
      return *m_src.getPixelPointer (x, y) < 128;
    }

  private:
    Pixels m_src;
  };

  /** Mask exclusion test functor.
  */
  struct BlackTest
  {
    explicit BlackTest (Pixels src) : m_src (src)
    {
    }

    inline bool operator () (int const x, int const y) const noexcept
    {
      return *m_src.getPixelPointer (x, y) > 127;
    }

  private:
    Pixels m_src;
  };

  /** Mask presence test functor.
  */
  struct AlphaTest
  {
    explicit AlphaTest (Pixels src) : m_src (src)
    {
    }

    inline bool operator () (int const x, int const y) const noexcept
    {
      return *m_src.getPixelPointer (x, y) != 0;
    }

  private:
    Pixels m_src;
  };

  //------------------------------------------------------------------------------

  /** Distance output to 8-bit unsigned.
  */
  struct OutputDistancePixels
  {
    OutputDistancePixels (Pixels dest, int radius)
      : m_dest (dest)
      , m_radius (radius)
      , m_radiusSquared (radius * radius)
    {
    }

    void operator () (int const x, int const y, double distance)
    {
      if (distance <= m_radiusSquared && distance > 0)
      {
        distance = sqrt (distance);

        *m_dest.getPixelPointer (x, y) = uint8 (255 * distance / m_radius + 0.5);
      }
      else
      {
        *m_dest.getPixelPointer (x, y) = 0;
      }
    }

  private:
    Pixels m_dest;
    int m_radius;
    int m_radiusSquared;
  };

  //------------------------------------------------------------------------------

  /** Distance output to a generic container.
  */
  template <class Map>
  struct OutputDistanceMap
  {
    typedef typename Map::Type Type;

    OutputDistanceMap (Map map, int radius)
      : m_map (map)
      , m_radius (radius)
      , m_radiusSquared (radius * radius)
    {
    }

    void operator () (int const x, int const y, double distance)
    {
      if (distance <= m_radiusSquared && distance > 0)
        m_map (x, y) = Type (std::sqrt (distance));
      else
        m_map (x, y) = 0;
    }

  private:
    Map m_map;
    int m_radius;
    int m_radiusSquared;
  };

  //------------------------------------------------------------------------------

  /** Distance output to a generic container.
  */
  template <class Map>
  struct OutputInverseDistanceMap
  {
    typedef typename Map::Type Type;

    OutputInverseDistanceMap (Map map, int radius)
      : m_map (map)
      , m_radius (radius)
      , m_radiusSquared (radius * radius)
    {
    }

    void operator () (int const x, int const y, double distance)
    {
      if (distance <= m_radiusSquared && distance > 0)
        m_map (x, y) = m_radius - Type (std::sqrt (distance));
      else
        m_map (x, y) = 0;
    }

  private:
    Map m_map;
    Type m_radius;
    Type m_radiusSquared;
  };

  //----------------------------------------------------------------------------
  // 
  // "A General Algorithm for Computing Distance Transforms in Linear Time"
  //  - A. Meijster, 2003
  //
  struct Meijster
  {
    struct EuclideanMetric
    {
      static inline int f (int x_i, int gi) noexcept
      {
        return (x_i*x_i)+gi*gi;
      }

      static inline int sep (int i, int u, int gi, int gu, int) noexcept
      {
        return (u*u - i*i + gu*gu - gi*gi) / (2*(u-i));
      }
    };

    struct ManhattanMetric
    {
      static inline int f (int x_i, int gi) noexcept
      {
        return abs (x_i) + gi;
      }

      static inline int sep (int i, int u, int gi, int gu, int inf) noexcept
      {
        if (gu >= gi + u - i)
          return inf;
        else if (gi > gu + u - i)
          return -inf;
        else
          return (gu - gi + u + i) / 2;
      }
    };

    struct ChessMetric
    {
      static inline int f (int x_i, int gi) noexcept
      {
        return jmax (abs (x_i), gi);
      }

      static inline int sep (int i, int u, int gi, int gu, int) noexcept
      {
        if (gi < gu)
          return jmax (i+gu, (i+u)/2);
        else
          return jmin (u-gi, (i+u)/2);
      }
    };

    template <class Functor, class BoolImage, class Metric>
    static void calculate (Functor f, BoolImage test, int const m, int const n, Metric metric)
    {
      std::vector <int> g (m * n);

      int const inf = m + n;

      // phase 1
      {
        for (int x = 0; x < m; ++x)
        {
          g [x] = test (x, 0) ? 0 : inf;

          // scan 1
          for (int y = 1; y < n; ++y)
          {
            int const ym = y*m;
            g [x+ym] = test (x, y) ? 0 : 1 + g [x+ym-m];
          }

          // scan 2
          for (int y = n-2; y >=0; --y)
          {
            int const ym = y*m;

            if (g [x+ym+m] < g [x+ym])
              g [x+ym] = 1 + g[x+ym+m];
          }
        }
      }

      // phase 2
      {
        std::vector <int> s (jmax (m, n));
        std::vector <int> t (jmax (m, n));

        for (int y = 0; y < n; ++y)
        {
          int q = 0;
          s [0] = 0;
          t [0] = 0;

          int const ym = y*m;

          // scan 3
          for (int u = 1; u < m; ++u)
          {
            while (q >= 0 && metric.f (t[q]-s[q], g[s[q]+ym]) > metric.f (t[q]-u, g[u+ym]))
              q--;

            if (q < 0)
            {
              q = 0;
              s [0] = u;
            }
            else
            {
              int const w = 1 + metric.sep (s[q], u, g[s[q]+ym], g[u+ym], inf);

              if (w < m)
              {
                ++q;
                s[q] = u;
                t[q] = w;
              }
            }
          }

          // scan 4
          for (int u = m-1; u >= 0; --u)
          {
            int const d = metric.f (u-s[q], g[s[q]+ym]) / 1;
            f (u, y, d);
            if (u == t[q])
              --q;
          }
        }
      }
    }
  };

  //----------------------------------------------------------------------------
  // 
  // "The 'Dead Recoking' Signed Distance Transform"
  // - G. Grevera, 2004
  //
  template <class T>
  static T distsq (T x, T y)
  {
    return sqrt (x*x + y*y);
  }

  struct Grevera
  {
    template <class Functor, class BoolImage, class Metric>
    static void calculate (Functor f, BoolImage test, int const X, int const Y, Metric)
    {
      typedef double T;
      typedef Point <int> P;

      T const inf = (T(X)*T(X) + T(Y)*T(Y));
      T const d1 = 1;
      T const d2 = sqrt (2.);

      Map2D <T> d (X, Y);
      Map2D <P> p (X, Y);

      for (int y = 0; y < Y; ++y)
      {
        for (int x = 0; x < X; ++x)
        {
          d (x, y) = inf;
          p (x, y) = P (0, 0);
        }
      }

      for (int y = 1; y < Y-1; ++y)
      {
        for (int x = 1; x < X-1; ++x)
        {
          bool const v = test (x, y);
          if ( (test (x-1,y) != v) ||
               (test (x+1,y) != v) ||
               (test (x,y-1) != v) ||
               (test (x,y+1) != v))
          {
            d (x, y) = 0;
            p (x, y) = P (x, y);
          }
        }
      }

      for (int y = 1; y < Y-2; ++y)
      {
        for (int x = 1; x < X-2; ++x)
        {
          if (d(x-1,y-1)+d2 < d(x,y))
          {
            p(x,y) = p(x-1,y-1);
            d(x,y) = distsq (T(x-p(x,y).getX()), T(y-p(x,y).getY()));
          }
          if (d(x,y-1)+d1 < d(x,y))
          {
            p(x,y) = p(x,y-1);
            d(x,y) = distsq (T(x-p(x,y).getX()), T(y-p(x,y).getY()));
          }
          if (d(x+1,y-1)+d2 < d(x,y))
          {
            p(x,y) = p(x+1,y-1);
            d(x,y) = distsq (T(x-p(x,y).getX()), T(y-p(x,y).getY()));
          }
          if (d(x-1,y)+d1 < d(x,y))
          {
            p(x,y) = p(x-1,y);
            d(x,y) = distsq (T(x-p(x,y).getX()), T(y-p(x,y).getY()));
          }
        }
      }

      for (int y = Y-2; y >= 1; --y)
      {
        for (int x = X-2; x >= 1; --x)
        {
          if (d(x+1,y)+d1 < d(x,y))
          {
            p(x,y) = p(x+1,y);
            d(x,y) = distsq (T(x-p(x,y).getX()), T(y-p(x,y).getY()));
          }
          if (d(x-1,y+1)+d2 < d(x,y))
          {
            p(x,y) = p(x-1,y+1);
            d(x,y) = distsq (T(x-p(x,y).getX()), T(y-p(x,y).getY()));
          }
          if (d (x,y+1)+d1 < d(x,y))
          {
            p(x,y) = p(x,y+1);
            d(x,y) = distsq (T(x-p(x,y).getX()), T(y-p(x,y).getY()));
          }
          if (d(x+1,y+1)+d2 < d(x,y))
          {
            p(x,y) = p(x+1,y+1);
            d(x,y) = distsq (T(x-p(x,y).getX()), T(y-p(x,y).getY()));
          }

          T dist = d(x,y);
          if (test (x, y))
            f (x, y, -dist * dist);
          else
            f (x, y, dist * dist);
        }
      }
    }
  };

  //----------------------------------------------------------------------------
  // 
  // "Efficient Euclidean Distance Transform Using Perpendicular Bisector Segmentation"
  // - J. Wang, Y. Tan, ?
  //
  struct WangTan
  {
    static int intersect (int ux, int vx, int du, int dv)
    {
      if (dv > du)
        return (dv - du) / (2 * (vx - ux));
      else
        return -2;
    }

    static int distance (int ux, int vx, int dv)
    {
      return ux * (ux - vx * 2) + dv;
    }

    template <class Functor, class BoolImage, class Metric>
    static void calculate (Functor f, BoolImage test, int const n, int const m, Metric)
    {
      int const inf = 1+n*n+m*m;
      Map2D <int> I (n, m);

      // stage 1
      for (int c = 0; c < n; ++c)
      {
        int mid = -1;
        int const cc = c*c;

        for (int r = 0; r < m; ++r)
        {
          if (test (c,r))
          {
            if (mid > -1)
              mid = (r + mid) /2;
            for (int r_ = mid + 1; r_ <= r; ++r_)
              I (c, r_) = (r_-r)*(r_-r) + cc;
            mid = r; // was x
          }
          else
          {
            if (mid == -1)
              I (c, r) = inf;
            else
              I (c, r) = (r-mid)*(r-mid) + cc;
          }
        }
      }

      // stage 2
      {
        std::vector <int> stack_c  (n);
        std::vector <int> stack_cx (n);
        std::vector <int> stack_g  (n);

        for (int r = 0; r < m; ++r)
        {
          int p = -1;

          for (int c = 0; c < n; ++c)
          {
            if (I (c, r) < inf)
            {
              for (;;)
              {
                int cx;

                if (p >= 0)
                {
                  cx = intersect (stack_c [p], c, stack_g [p], abs (I (c, r)));

                  if (cx == stack_cx [p])
                  {
                    --p;
                  }
                  else if (cx < stack_cx [p])
                  {
                    --p;
                    continue;
                  }
                  else if (cx >= (n - 1))
                  {
                    break;
                  }
                }
                else
                {
                  cx = -1;
                }

                ++p;
                stack_c  [p] = c;
                stack_cx [p] = cx;
                stack_g  [p] = I (c, r);
                break;
              }
            }
          }

          if (p < 0)
            return; // undefined!

          int c = 0;
          for (int k = 0; k <= p ; ++k)
          {
            int cx;
            if (k == p)
              cx = n-1;
            else
              cx = stack_cx [k+1];
            for (;c <= cx; ++c)
              I (c, r) = distance (c, stack_c [k], stack_g [k]);
          }
        }
      }

      // output
      for (int r = 0; r < m; ++r)
        for (int c = 0; c < n; ++c)
          f (c, r, I (c, r));
    }
  };
};

They are all extremely fast, and the first can be parallelized easily.

Any tips and pointers on accurately replicating Photoshop's "Layer Styles" would be most welcome.

[thevinn]

UPDATE!

We now know the exact algorithm:

1) Compute two values: "dilate" and "blur" measured in pixels. Both are dependent onthe values of both spread and size.

2) Apply "dilate" pixels of morphological dilation to the grayscale mask. The dilation needs to be a grayscale operator and not a binary operator. Photoshop does this by applying a Chamfer 5-7-11 distance transform using the technique exactly as described by G. Borgefors. The application of the dilation should produce fixed point integers with 8 bits of fractional precision (at no additional cost).

3) Apply "blur" pixels of box blur. The algorithm needs to support fractional radii. The easiest and fastest way to do this is to use four passes (rows, columns, rows, columns) with each pass transposing the output, with each pass at half the desired radius. The resulting intermediate image will be 8-bit fixed point.

4) Convert the intermediate image into a final alpha mask by taking each value and evaluating the 1-to-1 "contour" function. The output should be 8-bit unsigned integer in the range 0..255.

5) Fill the destination image with the drop shadow color, blend mode, and opacity using the final alpha mask.

This is all implemented in my cross platform demo application, which applies the drop shadow effect in real time as you change the controls:

https://github.com/vinniefalco/LayerEffects

[anthony]

In otherwords they use very specific ultra fast algorithms to go the job, where as IM is using very general multi-purpose algorithms.
The results however should be almost identical, with with a factional radial dialtion method basied on a distance algorithm.

NOTE the distance algorithm built into IM is fractional and uses a fast 2 pass method (though still uses a 2-D kernel).

Hmmm the names are different. I named the mtrics for the distance metric while the papers is naming then for the algorithm designer
Chessboard is Chebyshev
CityBlock is Manhattan
And Montanari is Euclidean (knight)

Chamfor appears to be a special integer form of radius 2 Eulcidean using a distance kernel of the form
- 11 - 11 -
11 7 5 7 11
- 5 0 5 -
11 7 5 7 11
- 11 - 11 -

NOTE normal distance kernels in IM use 100 as the orgongonal neighbour distance scale, where chamfor uses a scale of '5'
It may be a interesting 'specialized' distance kernel appropriate to images needing HIGH distances, and a good interger distance results.

Okay.. i rambling.. I'll be quiet now.

Hmmm I have added some notes in IM examples about the Chamfor kernel (whcih is really a set of kernels) but it is not implemented.
http://www.imagemagick.org/Usage/morphology/#chamfer