Skip to main content

plotting - Custom interpolation algorithm


Suppose I have a list to create a list density plot, the list format is:


list = {{0.5, 0.5, 0.0}, {-0.5, 0.5, 1.0}, {1.0, -0.45, 0.0}, {0.0, 0.0, 0.0}}

The size of the list is about 10^5 elements


The plot will show a point or nothing at all. The real data produces this picture:


enter image description here


My goal is to make solid lines so that the shapes look smooth, however, I also need a list containing all the points that make the smoothed out picture. Since I know the distance of the mesh, I was thinking about a procedure that would go through the whole grid searching for neighboring points and if there exists a point with {x, y, 1} would fill in the blanks. I'm stuck in the criterion to use for this though. Any thoughts?




Answer



I imported your posted data:


raw = Import["array1.CSV"];

Then selected one of the features, and reduced it to two-dimensional, since the third dimension is uniformly equal to $1$ for the points constituting the feature:


single = Cases[raw, {a_, b_, 1.} :> {a, b} /; (-0.2 < a < 0.2 && 0.2 < b)];

ListPlot[
single,
PlotRange -> {0.2, 0.5}, AspectRatio -> Automatic

]

plot


I isolate the oval feature manually:


oval = Cases[single, {a_, b_} /; (-0.07 < a < 0.07 && 0.35 < b < 0.472)];

ListPlot[
oval,
PlotRange -> {0.39, 0.48}, AspectRatio -> Automatic
]


oval points


Then I fit an axis-aligned ellipse to it, by generating an appropriate squared distance function to be minimized, then making the assumption that the abscissa of the center is $0$, which seems reasonable given the symmetry of the system, and imposing reasonable constraints on the other parameters obtained from inspection of the graph:


obj[a_, b_, xc_, yc_] = Total[
Simplify[
SquaredEuclideanDistance[{#1, #2}, (#1 - xc)^2/a^2 + (#2 - yc)^2/b^2 - 1] & @@@ oval,
_ ∈ Reals
]
];


minPars = FindMinimum[{obj[a, b, 0, yc], 0.42 < yc < 0.46}, {{a, 0.5}, {b, 0.5}, {yc, 0.44}}]

(* {14.2803, {a -> 0.0518811, b -> 0.0312114, yc -> 0.439966}} *)

It's a pretty good fit, despite the few stray points we had to tolerate:


cplot = ContourPlot[
Evaluate[((x - 0)^2/a^2 + (y - yc)^2/b^2 == 1) /. First@Rest@minPars],
{x, -0.08, 0.08}, {y, 0.4, 0.48},
Epilog -> Point[oval], AspectRatio -> Automatic
]


fit oval and points


Now, we can extract the line shape from the contour plot results:


ovalLine = First@Cases[Normal@cplot, _Line, Infinity];



I then carefully select the points at the periphery of your feature, through somewhat laborious manual filtering:


rdf = RegionDistance[ovalLine];

externalPoints = Join[

DeleteCases[single, pt_ /; rdf[pt] < 0.00715],
Complement[
Select[single, #[[2]] > 0.4713 &],
MinimalBy[Select[single, #[[2]] > 0.4713 &], Abs@#[[1]] &, 2]
]
];

Graphics@Point@externalPoints

external



I then use the fantastic alphaShapes2DC function proposed by RunnyKine in this answer, to generate a "concave hull" of those points:


alphaShapes2DC[externalPoints, .10]

alpha shape


In fact, I will turn the 2D region returned by that function into a boundary mesh region, with some formatting for appearances only:


reg = BoundaryDiscretizeRegion[
alphaShapes2DC[externalPoints, .10],
MeshCellStyle -> {{1, All} -> Directive[Thick, Red]},
PlotTheme -> "Lines"
];


formatted alpha shape




Let's now put it all together:


Show[reg, Graphics[{Thick, Red, ovalLine}]]

both lines


And here's a comparison to the original points:


Show[reg, Graphics[{Thick, Red, ovalLine, Black, Point[single]}]]


with points




And finally let's generate polygons and rotate them around the axis to generate the final shape:


allShapes = NestList[
GeometricTransformation[#, RotationTransform[2 Pi/3, {0, 0}]] &,
MeshPrimitives[reg, 2],
2
];

allOvals = NestList[

GeometricTransformation[#, RotationTransform[2 Pi/3, {0, 0}]] &,
Polygon @@ ovalLine,
2
];

Graphics[{FaceForm[None], EdgeForm[{Thick, Red}], allOvals, allShapes}]

final


Comments

Popular posts from this blog

plotting - How to draw lines between specified dots on ListPlot?

I would like to create a plot where I have unconnected dots and some connected. So far, I have figured out how to draw the dots. My code is the following: ListPlot[{{1, 1}, {2, 2}, {3, 3}, {4, 4}, {1, 4}, {2, 5}, {3, 6}, {4, 7}, {1, 7}, {2, 8}, {3, 9}, {4, 10}, {1, 10}, {2, 11}, {3, 12}, {4,13}, {2.5, 7}}, Ticks -> {{1, 2, 3, 4}, None}, AxesStyle -> Thin, TicksStyle -> Directive[Black, Bold, 12], Mesh -> Full] I have thought using ListLinePlot command, but I don't know how to specify to the command to draw only selected lines between the dots. Do have any suggestions/hints on how to do that? Thank you. Answer One possibility would be to use Epilog with Line : ListPlot[ {{1, 1}, {2, 2}, {3, 3}, {4, 4}, {1, 4}, {2, 5}, {3, 6}, {4, 7}, {1, 7}, {2, 8}, {3, 9}, {4, 10}, {1, 10}, {2, 11}, {3, 12}, {4, 13}, {2.5, 7}}, Ticks -> {{1, 2, 3, 4}, None}, AxesStyle -> Thin, TicksStyle -> Directive[Black, Bold, 12], Mesh -> Full, Epilog -> { Line[ ...

calculus and analysis - Integral of the Sinc product

Let us consider the following integral $$ B_n = \int_0^\infty \prod_{k=1,3,5,\dots}^n\frac{\sin (x/k)}{x/k}dx $$ By definition, Sinc[x] == Sin[x]/x , therefore B[n_?OddQ] := Integrate[Product[Sinc[x/k], {k, 1, n, 2}], {x, 0, Infinity}] B /@ Range[1, 13, 2] {Pi/2, Pi/2, Pi/2, Pi/2, Pi/2, Pi/2, Pi/2} OK, everything is fine. But... B[15] 467807924713440738696537864469 Pi/935615849440640907310521750000 What's going on? Answer As Eckhard wrote in comments B[n] is the n-th Borwein integral . (The letter B was not accidental :) ) This funny properties of Borwein integrals is related to the Fourier transform of Sinc function FourierTransform[Sinc[x], x, k] 1/2 Sqrt[Pi/2] (Sign[1 - k] + Sign[1 + k]) Plot[%, {k, -2, 2}, Filling -> 0] which is the box function. The result is $\pi/2$ while the sum $$ 1/3+1/5+\dots+1/n $$ B_n = \frac{\pi}{2} - \pi \bigg(\sum_{k=3,5,\ldots}^n\frac{1}{k} -1\biggr)^\frac{n-1}{2}\prod_{k=3,5,\dots}^n\frac{k}{k-1}. $$ For $n=15$ it is equal to $$ \frac{46780...

list manipulation - Selecting multiple columns from a matrix?

Sample data: data = { {{2013, 1, 1}, 24.13, 167.67, 231.82}, {{2013, 1, 2}, 32.15, 170.92, 225.99}, {{2013, 1, 3}, 35.43, 172.68, 221.67}, {{2013, 1, 4}, 36.73, 173.05, 218.32}, {{2013, 1, 5}, 58.19, 165.96, 197.05}, {{2013, 1, 6}, 69.99, 163.50, 187.52}, {{2013, 1, 7}, 71.37, 154.21, 175.58}, {{2013, 1, 8}, 72.51, 149.66, 163.25}}; I want a DateListPlot with three graphs, so for a matrix formed by columns 1 and 2, one for columns 1 and 3, and 1 for columns 1 and 4. At the moment I'm using this code: data2 = Transpose[{data[[All, 1]], data[[All, 2]]}]; data3 = Transpose[{data[[All, 1]], data[[All, 3]]}]; data4 = Transpose[{data[[All, 1]], data[[All, 4]]}]; DateListPlot[{data2, data3, data4}, Joined -> True, Filling -> {3 -> {1}}] but I have a hunch that this can be done more efficiently. I don't like the Transpose s in particular. Any ideas? edit (for extra credit) What if I need to multiply the second column by 2, which in my solution is simp...