# Tag Info

18

EDIT: graphical alignment issue corrected. I think this has promise. The centering is almost perfect, and while the code may not be well tuned or optimized it is quite fast: 0.008736 seconds on my machine. It works by attempting to find the center of each white "blob" and then averaging those positions. img = Import["http://i.stack.imgur.com/i050B.png"]; ...

18

Your image: img = Import["http://i.stack.imgur.com/i050B.png"]; Here is the circle's radius and center: crcl = ComponentMeasurements[RegionBinarize[Dilation[img, 3], {{320, 240}}, 0.3], {"Centroid", "EquivalentDiskRadius"}]; // AbsoluteTiming crcl {0.126007, Null} {1 -> {{284.448, 241.873}, 109.256}} And here how precise it is: Show[img, ...

16

$Version (* "8.0 for Microsoft Windows (64-bit) (October 7, 2011)" *) Direct attack fails: Timing[Convolve[Sinc[x], Exp[-x^2], x, y]] (* Out= {59.296, Convolve[Sinc[x], E^-x^2, x, y]} *) or, equivalently, Timing[Integrate[Sinc[x] Exp[-(x - y)^2], {x, -∞, ∞}] ]$\left\{49.92,\int_{-\infty }^{\infty } e^{-(x-y)^2} \text{Sinc}[x] \, dx\right\}$... 15 All the image samples you posted are basically lines pointing towards a common center. We can turn that into a simple mathematical model: Let's say every gradient in the image is normal to a line pointing towards the center. Find the center with the least squared error. So the error term would be: squaredError = ({cx - x, cy - y}.{gx, gy})^2; where cx/cy ... 14 I hope I see the essence here. You are interested in the convolution of an interpolated function with a Gauss function Your underlying data has regular spacings in x-direction and the convolution with a Gaussian is extremely fast implemented in GaussianFilter for discrete data. Why are you making it so complicated when the only thing you have to do is to ... 14 Here I implemented three different versions of Gaussian filtering (for periodic data). It took me a while to adjust the constants and still some of them might be wrong. Prepare the Gaussian kernel n = 200000; σ = .1; t = Subdivide[-1. Pi, 1. Pi, n - 1]; ker = 1/Sqrt[2 Pi]/ σ Exp[-(t/σ)^2/2]; ker = Join[ker[[Quotient[n,2] + 1 ;;]], ker[[;; Quotient[n,2]]]]; ... 13 In line with the OPs request for a comparison of several different methods here's a comparison of six different ways to filter (or convolve) a data set x with a kernel h: convolution, correlation, the frequency domain method, a direct time-domain method such as might be programmed in C or Java, and a vectorized version such as would be common in Matlab or ... 12 The functions do not have a finite area, so they cannot be real distributions as your title claims they are. Let's change them a bit so they have area 1. f[x_] = (1/k) Exp[-x/k] UnitStep[x]; g[x_] = (1/p) Exp[-x/p] UnitStep[x]; Integrate[f[x], {x, -∞, ∞}] ConditionalExpression[1, Re[1/k] > 0] The convolution: Convolve[f[x], g[x], x, y] which ... 12 As mentioned in a comment by ciao, blind deconvolution is no easy task. However, here are some pointers for estimating the kernel parameters. These can equally be applied to denoising an image rather than deblurring. I am assuming that the form of the kernel is known here, rather than blind deconvolution (Wikipedia) methods such as Maximum A Posteriori (... 11 While it's tempting to attribute the errors you're observing to floating-point errors due to zeros in the DFT of the window, this is actually not the case here: Window[width_, x_] := UnitStep[x + width/2] UnitStep[width/2 - x]; Test[x_] := UnitStep[x] UnitStep[count - x] Sin[1/10 x]^2; tempWindow = Table[Window[20, x], {x, -count/2, count/2}]; {Min[Abs[... 10 In this case, the following works just fine: LeastSquares[ Transpose@Through[{RotateLeft, Identity, RotateRight}[data1]], conv ] {0.333333, 0.333333, 0.333333} More general with kernel length 2 k + 1: a = RandomReal[{-1, 1}, 100]; k = 20; ker = RandomReal[{-1, 1}, 2 k + 1]; b = ListConvolve[ker, a, k + 1]; ker1 = LeastSquares[ Transpose[... 9 Unfortunately, MXNet does not support 3D convolutions with dilations yet. This can be seen in the MXNet source for convolution. 9 Consider the following code f[n_, p_] := n^p g[n_, p_] := n*p DirichletConvolve[f[n, p], g[n, p], n, 4] First, we define two functions f and g. Then we compute their Dirichlet convolution. The third argument in the Dirichlet convolution tells us that n is the function argument for which we want to do the convolution. p on the other hand is a parameter ... 8 I think you are trying to do the following: Clear[h]; h[n_] /; n >= 1 := -((I^-n (-1 + I^n)^2)/(n^2 Pi^2)) h = 1/4; Clear[y]; y[m_] := DiscreteConvolve[h[n], DiscreteDelta[n], n, m] y (* ==> 1/4 *) Here corrected your definition of y so it uses m as the function argument, but the important part is to understand why even with this ... 8 Here is a way to solve this problem using the convolution theorem: l = Assuming[{γ > 0 && σ > 0 && μ > 0 && k ∈ Reals}, FourierTransform[PDF[CauchyDistribution[μ, γ], x], x, k] ]$\frac{\left(\theta (-k) e^{2 \gamma k}+\theta (k)\right) e^{-k (\gamma -i \mu )}}{\sqrt{2 \pi }}$g = Assuming[{γ > 0 &&... 8 That does not explain why your solution does not work (I'm pretty sure it is because of your definition of stripe), but if you replace stripe with: stripe[x_, d_] := HeavisidePi[Mod[x/(2 d), 1]] you get: conv[x_] = Convolve[UnitStep[tau + 5] stripe[tau, 1] UnitStep[5 - tau], gaussian[tau, 0, 0.5], tau, x] (* -0.5 (1. Erf[2.82843 - 1.41421 x] + 1. ... 7 I'm not sure you are handling the top and left edges in the way you really want; they work with the second rather than first elements being treated as "middle" twice, with first elements not treated that way at all. Here is code that does not do that, hence gives different results than yours on top and left edges. It is around two orders of magnitude faster.... 7 Related to your question: If you wish to do general cross-correlation (not with a fixed lag of 1), you can use the following functions: xcorrlag[dat1_, dat2_] := Position[#, Max@#] &[xcorr[dat1, dat2]] - Length[dat1] xcorr[dat1_, dat2_] := ListCorrelate[dat1, dat2, {-1, 1}, 0] xcorr[dat_] := xcorr[dat, dat] Here, xcorr called with a single Listas an ... 7 We consider h = Exp[g]. The main idea is to calculate the Fourier transform of h, expand it into a series in powers of v, and then transform back Here we go. The normal distribution is explicitly f = 1/(v Sqrt[2 \[Pi]]) Exp[-(x - z)^2/(2 v^2)]; Now the Fourier transformation of h acts only on f, with the result ft = FourierTransform[f, x, t, ... 7 I wasn't sure how to set it up with the standard fit functions, so I just rolled my own least squares. (It probably can be done, but I had an inkling I might want to have greater control over the computation. I hope it helps.) I start by defining m1 and m2 on your data: Clear[m1, m2]; (m1[t_] /; t == First@# = Last@#) & /@ m1data; (m2[t_] /; t == ... 7 It seems, that MMA can't do that integral analytically (working with Version 8.0). Also Convolvedidn't do the job. Do it numerically. By the way, don't use v', because it is interpreted as Derivative and regard, the correct definition is gl[v - vs, v0, Lw], not vs-v. gd[v_, v0_, Dw_] = Sqrt[(4*Log/Pi)]*(1/Dw)*Exp[-4*Log*((v - v0)/Dw)^2]; gl[v_, ... 7 Let's focus on the part that is giving you trouble: the term$\sum_{k=0}^{m-1}k^t(m-k)^t\$, which is written in code as: Sum[k^t (m - k)^t, {k, 0, m - 1}] For example, with m=30, you get: Your stated goal is to express this in term of a discrete convolution. Here's one way to do that. Let kt = Range[0, 30]^t; be one of the signals. Then the sum is ...

7

It's hard to know for sure, but one way to test for caching is to apply a single command to lots of data sets, or to apply the command separately to each set. For instance: n = 5000; data = RandomReal[{-1, 1}, {n, 10000}]; GaussianFilter[#, 100] & /@ data; // AbsoluteTiming Do[GaussianFilter[data[[i]], 100], {i, n}] // AbsoluteTiming Do[GaussianFilter[...

7

Not sure if I've guessed definition of conv2 accurately, but the following does reproduce the output of MATLAB. (Tested in Octave. ) conv2[A_?NumericQ, B_] := conv2[{{A}}, B] conv2[A_?VectorQ, B_] := conv2[{A}, B] conv2[A_, B_] := ListConvolve[A, B, {1, -1}, 0] filter = 1; F = conv2[{1, 2, 1}, List /@ {1, 2, 1}]/16.; Nest[conv2[#, F] &, filter, 3] For ...

7

(*Defining Tent Funtion*) Tri[x_, wT_] := Piecewise[{ {(1/wT^2)*x + (1/wT), 0 > x > -wT}, {-(1/wT^2)*x + (1/wT), 0 < x < wT}}]; (*Defining Top Hat Funtion*) Hat[x_, wH_] := Piecewise[{ {(1/(2*wH)), 0 > x > -wH}, {(1/(2*wH)), 0 < x < wH}}]; (*Plotting& Manipulate Function*) Manipulate[ s1 = Switch[fx, "Tent ...

6

You can do the cross-correlation between two sequences using ListCorrelate a = {1, 2, 3, 4, 5}; b = {5, 4, 3, 2, 1}; ListCorrelate[a, b, 1] This does circular correlation, so you may want to look at the options to get the exact calculation you are looking for.

6

(1) Convolution with Plus and Times can be done via FFT. (2) The overall speed complexity cannot be less than the size-of-result complexity. Point (1) might help to explain why the standard ListConvolve is fast under most circumstances. Point (2) on the examples in this post should help to explain why LC2 is likely to be slow. To make this clear we can ...

6

Automatic is not the same as "Varying". Example 1 This works: conv1 = NetInitialize@ConvolutionLayer[3, 2, "Input" -> {3, Automatic, Automatic}]; conv1@NetEncoder[{"Image", {128, 128}}]@RandomImage[1, {128, 128}, ColorSpace -> "RGB"] This doesn't work at the stage of layer creation: conv2 = NetInitialize@ConvolutionLayer[3, 2, "Input" -> {3, "...

6

First, you need to define the convolution of the two functions. Then plot it: con[y_] := Convolve[HeavisidePi[x], HeavisidePi[x], x, y]; Manipulate[Show[Plot[HeavisidePi[y], {y, -3, 3}, Exclusions -> None, Filling -> Bottom], Plot[HeavisidePi[x - t], {x, -3, 3}, Exclusions -> None, PlotStyle -> Green, Filling -> Bottom], Plot[...

5

Sparse multiplication with a circular matrix corresponds to a convolution; on a trivial example let us compare: matrix = SparseArray[{Band[{1, 1}] -> 2, Band[{1, 2}] -> 1, Band[{2, 1}] -> 1}, {15, 15}]; vec = SparseArray[5 -> x, 15]; matrix.vec // Normal (* ==> {0, 0, 0, x, 2 x, x, 0, 0, 0, 0, 0, 0, 0, 0, 0} *) versus a = SparseArray[5 -&...

Only top voted, non community-wiki answers of a minimum length are eligible