1
$\begingroup$

I'm currently finishing some research from a internship I did last summer and part of the project involves calculating the centroids of coronal holes from synoptic maps. If you aren't aware of what these are, they are basically holes or regions of the corona, an upper layer of the suns atmosphere, where it is less dense. The purpose of this is so I can, track how the motion of these coronal holes change over time. The way this is done is by taking a certain latitude band, say 10 degrees south to 10 degrees north, finding the centroids for coronal holes after each carrington rotation, solar month, and finding the slope with all of these points. In general, calculating these centroids is not an easy task because in general the shapes of these regions, on solar maps, are non uniform and can take on any flat 2D shape. However last summer I did find a solution that should work, at least theoretically. Essentially my idea was to generate parametric plots that look identical to the images of each coronal hole and then use green's theorem, with the parametric equations, to calculate the centroids through line integration. Calculus 3 was still fresh in my mind at that point and I thought, if we can use Green's Theorem to calculate area, why not centroids as well? After doing some research I found that these line integrals do in fact exist and I went through the derivations to prove to myself that the line integrals were equivalent to the double integrals used to calculate centroids. If you are wondering what these line integrals are, here they are.

enter image description here

And here are the double integrals that the line integrals are equivalent to.

enter image description here

My mentors had another method that they had planned for me to use on IDL but encouraged me to try applying my method instead. I ended up working on Mathematica instead, because of useful code I found on other forums. From there I added some of my own code and did a bunch of tinkering over the course of several weeks. At the time, I did not know if this was the most efficient method computationally and I still don't but I would still like to get this to be fully functional. I got the method to work for the most part, but at times I'm still having issues getting an accurate answer. Sometimes it gives an answer that looks reasonable but at other times it gives an answer that is obviously completely off just by eyeballing it. I believe that the issue solely lies in convergence issues with the integrals but I will show all of my code, just in case some of the issues lie elsewhere.

param[x_, m_, t_] := 
 Module[{f, n = Length[x], nf}, 
  f = Chop[Fourier[x]][[;; Ceiling[Length[x]/2]]];
  nf = Length[f];
  Total[Rationalize[
     2 Abs[f]/Sqrt[n] Sin[
       Pi/2 - Arg[f] + 2. Pi Range[0, nf - 1] t], .01][[;; 
     Min[m, nf]]]]]

tocurve[Line[data_], m_, t_] := param[#, m, t] & /@ Transpose[data]
img = \!\(\*
GraphicsBox[
TagBox[RasterBox[CompressedData["
1:eJzt229sE3Ucx/HP/f/dXa/XrmUmiFl8hCgzMWbBTDfnnH+GmyYVE+IziJMY
khG3ZUYhZiISM0WmENQwUNRoGNncEP8skoFREQMD94dtOjIyKZtUxv5gO2XZ
6ukjn5ioafzu4PtK2vTaJ+82l+v3d71ev3pd7HEZQLXw7mKrnrqrqmrV0w+H
vI1HKqufWFNZ8VhpZU3Fmoqq21Yr3pMLvds6CfjjcZoxxhhjjDHGGGOMZU7s
7r22uhvYodhNktMCvC9rb1BHMZ95fdu3QIy6ghHLW9oHHIuGxywxCHwXDl2S
EAe+iYpG6jTmG0CdJn9RU9tFHZJh5eWtwHvAx0Bjff2V9u4y697ir71DRyR8
2XWms7LStjkDpBw97Wjp667pp65jvrFs2X6B4fz8buqQDAPeCVs/Rs1kUBsU
yn7qnHntjqI24CjQD3QBA8CIrkxYxoQuex/j996ShzqQ+QbQrWD4xps/pw7J
mA11nd5Mku2OuWIibMRltFIXzXfPbj+7qaF/8yvx+q2Jl172FjhdhjYOTKr2
jOWe9Y4nu976gbqRzXfLlw/YYkRYvwAnnnvhDHVOxtx6i7fA+Sxk9oQCHdQt
vuQEhlR9XDaSemhWUfoltfW1LbxaZH+r9J7jQKfrJAN2WlHGluZ+SV2UYSWF
zcBK6gpfAj6QtQSQgpGG9psuj0jYRx3F5ilghYrjAZECzgtz3HCTwMmcRV33
l3QW33mw6PZ26kBGTNL6nKw5SYlbkcvAsdwlV85CmGUc0GaK85HorOmmJHVY
038ylOSfZ+EOALtzb/h08WK+2ODqJUmfCOecKk4BTQ/EDlHnsPlOlt82xQVZ
vhhwx5yslKqe1kV79qIPc3JaVGwPiUbXraFuZAQCzpshpztg93vfLPmF/IMO
+0eAoVBkzs2aAno2bkr89aWCvJ1UVYxQMNhqqwlb8naJ9kdjR6hzmD/UbRyV
5IShz2rWpCQdamgYpC5ixCAaVNEbDU5HrZGVZT3UOcw3gL1BZ9wUSVkZta2p
J9f+TF3EyNiRFm8UkcWZcE4KaJWkHdRFzDc2VJ821N6g8GaSSU2ZjthzQNuW
zX3UXez/lpfXBOwzjSHLuRRwR4Fda2tPUEcxPwGaLT1uKNOqkpakX01lqKSw
kzqKESgu+8owTmrKcLY1xtfPs//gmfUXJPmiok6ZYiYUnPb2ond38lhyNSoo
/QjoEPKAqw9u28p/7mP/WvmDRxRlUNVGDRFfsOAcsIe6iNEoKtujKs3Xhvts
6TB1C/OfJTc9DzTZdqdu9GrGKdM8KsuvUkcxGgX3vQjURmxvPmmgbmE+E134
EBATdrOwDgMtmt7hugeB9dRdjDHGGGOMMcYYY4xl2O/DK24g
"], {{0, 26}, {367, 0}}, {0, 255},
ColorFunction->RGBColor],
BoxForm`ImageTag["Byte", ColorSpace -> "RGB", Interleaving -> True],
Selectable->False],
DefaultBaseStyle->"ImageGraphics",
ImageSizeRaw->{367, 26},
PlotRange->{{0, 367}, {0, 26}}]\)
img = Binarize[img~ColorConvert~"Grayscale"];
lines = Cases[
   Normal@ListContourPlot[Reverse@ImageData[img], 
     Contours -> {0.5}], _Line, -1];

ParametricPlot[
 Evaluate[tocurve[#, 100000000, t] & /@ lines], {t, 0, 1}, 
 Frame -> True, Axes -> False, ImageSize -> Large]
ParametricEquations = Evaluate[tocurve[#, 100000000, t] & /@ lines]
n = Length[ParametricEquations]
img2 = Binarize[img~ColorConvert~"Grayscale"];

bsFs = Cases[
   Normal@ListContourPlot[Reverse@ImageData[img2], Contours -> {0.5}],
    Line[x_] :> BSplineFunction[x], All];

ParametricPlot[Evaluate[Through@bsFs@t], {t, 0, 1}, Frame -> True, 
 Axes -> False, ImageSize -> Large, Axes -> false, Frame -> True]
stepsize = 0.0001
list = Table[
  N[Evaluate[Through@bsFs@t] /. t -> i], {t, 0, 1, stepsize}]
xmatrix = 
 Table[Table[list[[i, j, 1]], {i, 1, (1/stepsize) + 1}], {j, 1, n}]
ymatrix = 
 Table[Table[list[[i, j, 2]], {i, 1, (1/stepsize) + 1}], {j, 1, n}]
Minx = Table[xmatrix[[i, Ordering[xmatrix[[i]], 1]]], {i, 1, n}]
Miny = Table[ymatrix[[i, Ordering[ymatrix[[i]], 1]]], {i, 1, n}]
x = Table[ParametricEquations[[i, 1]], {i, 1, n}]
y = Table[ParametricEquations[[i, 2]], {i, 1, n}]
Mx = Table[NMinValue[{x[[i]], 0 <= t <= 1}, t], {i, 1, n}]
My = Table[NMinValue[{y[[i]], 0 <= t <= 1}, t], {i, 1, n}]
NewParametricEquations = 
 Table[{x[[i]] - (Mx[[i]] - Minx[[i]]), 
   y[[i]] - (My[[i]] - Miny[[i]])}, {i, 1, n}]
xCentroids = 
 Table[NIntegrate[
    NewParametricEquations[[i, 1]]^2*
     D[NewParametricEquations[[i, 2]], t], {t, 0, 1}]/
   NIntegrate[
    2*NewParametricEquations[[i, 1]]*
     D[NewParametricEquations[[i, 2]], t], {t, 0, 1}], {i, 1, n}]
yCentroids = 
 Table[NIntegrate[-NewParametricEquations[[i, 2]]^2*
     D[NewParametricEquations[[i, 1]], t], {t, 0, 1}]/
   NIntegrate[-2*NewParametricEquations[[i, 2]]*
     D[NewParametricEquations[[i, 1]], t], {t, 0, 1}], {i, 1, n}]
Centroids = Table[{xCentroids[[i]], yCentroids[[i]]}, {i, 1, n}]
plot = ParametricPlot[NewParametricEquations, {t, 0, 1}, 
   ImageSize -> Large, Axes -> False, Frame -> True, PlotRange -> All];
centroidpoints = 
  ListPlot[Table[{First[xCentroids[[i]]], First[yCentroids[[i]]]}, {i,
       1, n}] -> Range[n], LabelingFunction -> Bottom, 
   PlotStyle -> {Red, PointSize[Medium]}];
coronalholepoints = 
  ListPlot[Table[{First[N[NewParametricEquations[[i, 1]] /. t -> 0]], 
      First[N[NewParametricEquations[[i, 2]] /. t -> 0]]}, {i, 1, 
      n}] -> Range[n], LabelingFunction -> Left, 
   PlotStyle -> {Blue, PointSize[Medium]}];
Show[plot, centroidpoints, coronalholepoints]

Here is an example of an image with some coronal holes I attempted to calculate the centroids of. All other colors have been filtered out in order to generate accurate parametric equations for each coronal hole (the blue regions in this case, which indicates that the coronal holes have a positive polarity). enter image description here

And here is an image of the parametric plot of all the coronal holes, along with red dots which represent the calculated centroids. enter image description here As you can see, the first two centroids look reasonable but the last one is way off and isn't even inside the coronal hole. As a reference, here is a parametric equations that Mathematica generated for the coronal holes.

{{{24.7056 - 1/45 Sin[1/11 - (3349 t)/13] - 
    1/77 Sin[3/8 - (754 t)/3] - 1/23 Sin[8/9 - (3217 t)/16] - 
    1/76 Sin[1/18 - (377 t)/2] - 1/28 Sin[1/4 - (1640 t)/9] - 
    1/30 Sin[7/13 - (1797 t)/11] - 1/53 Sin[2/9 - (1885 t)/12] - 
    1/23 Sin[22/15 - (4191 t)/29] - 1/28 Sin[5/9 - (2111 t)/16] - 
    1/22 Sin[7/5 - (2199 t)/25] - 1/16 Sin[5/4 - (1552 t)/19] - 
    1/9 Sin[1/4 - (553 t)/8] - 1/11 Sin[8/9 - (377 t)/6] - 
    1/8 Sin[13/10 - (553 t)/11] - 59/15 Sin[13/12 - (88 t)/7] + 
    1/5 Sin[(1891 t)/43] + 19/5 Sin[17/5 + (44 t)/7] + 
    5/4 Sin[103/26 + (132 t)/7] + 2/3 Sin[25/14 + (201 t)/8] + 
    6/13 Sin[47/14 + (377 t)/12] + 1/7 Sin[13/4 + (377 t)/10] + 
    1/9 Sin[25/6 + (509 t)/9] + 1/7 Sin[5/14 + (377 t)/5] + 
    1/24 Sin[23/14 + (377 t)/4] + 1/24 Sin[45/13 + (1307 t)/13] + 
    1/16 Sin[8/5 + (1175 t)/11] + 1/11 Sin[94/21 + (1131 t)/10] + 
    1/71 Sin[31/7 + (955 t)/8] + 1/25 Sin[13/5 + (377 t)/3] + 
    1/27 Sin[3/14 + (1244 t)/9] + 1/19 Sin[5/4 + (754 t)/5] + 
    1/33 Sin[65/22 + (2287 t)/13] + 1/23 Sin[58/19 + (1753 t)/9] + 
    1/27 Sin[3/2 + (3525 t)/17] + 1/55 Sin[20/7 + (1709 t)/8] + 
    1/65 Sin[14/3 + (2419 t)/11] + 1/28 Sin[10/19 + (1131 t)/5] + 
    1/85 Sin[49/11 + (3952 t)/17] + 1/43 Sin[10/3 + (3104 t)/13] + 
    1/33 Sin[5/11 + (5391 t)/22] + 1/56 Sin[7/2 + (2375 t)/9] + 
    1/43 Sin[14/15 + (3041 t)/11] + 1/92 Sin[41/11 + (1131 t)/4] + 
    1/80 Sin[32/9 + (2953 t)/10]}, {14.8678 - 
    1/86 Sin[5/4 - (1885 t)/6] - 1/83 Sin[16/15 - (3349 t)/13] - 
    1/21 Sin[6/5 - (3217 t)/16] - 1/21 Sin[3/5 - (1640 t)/9] - 
    1/46 Sin[7/11 - (1797 t)/11] - 1/36 Sin[7/8 - (1244 t)/9] - 
    1/11 Sin[4/7 - (2111 t)/16] - 1/44 Sin[17/18 - (955 t)/8] - 
    1/34 Sin[7/6 - (1307 t)/13] - 1/15 Sin[1 - (2199 t)/25] - 
    1/7 Sin[23/15 - (553 t)/11] - 3/13 Sin[23/15 - (377 t)/10] - 
    2/5 Sin[7/10 - (377 t)/12] - 8/11 Sin[4/3 - (132 t)/7] - 
    10/7 Sin[5/6 - (88 t)/7] - 145/18 Sin[13/9 - (44 t)/7] + 
    4/11 Sin[12/11 + (201 t)/8] + 2/9 Sin[1/3 + (1627 t)/37] + 
    1/8 Sin[14/5 + (509 t)/9] + 1/24 Sin[137/46 + (377 t)/6] + 
    1/10 Sin[35/12 + (553 t)/8] + 1/66 Sin[5/8 + (377 t)/5] + 
    1/9 Sin[57/14 + (1307 t)/16] + 1/30 Sin[5/6 + (1175 t)/11] + 
    1/22 Sin[50/11 + (1131 t)/10] + 1/21 Sin[8/3 + (377 t)/3] + 
    1/30 Sin[22/5 + (3324 t)/23] + 1/22 Sin[26/11 + (754 t)/5] + 
    1/27 Sin[34/15 + (1885 t)/12] + 1/24 Sin[19/7 + (1866 t)/11] + 
    1/69 Sin[23/8 + (2287 t)/13] + 1/23 Sin[1/5 + (377 t)/2] + 
    1/30 Sin[11/3 + (1753 t)/9] + 1/21 Sin[37/11 + (1709 t)/8] + 
    1/64 Sin[29/7 + (2419 t)/11] + 1/38 Sin[10/19 + (1131 t)/5] + 
    1/31 Sin[4/5 + (5391 t)/22] + 1/36 Sin[29/8 + (2375 t)/9] + 
    1/86 Sin[5/7 + (3041 t)/11] + 
    1/64 Sin[48/11 + (1131 t)/4]}}, {{217.089 - 
    1/41 Sin[36/35 - (2375 t)/9] - 1/65 Sin[2/3 - (3952 t)/17] - 
    1/17 Sin[1/30 - (3217 t)/16] - 1/18 Sin[31/21 - (2287 t)/13] - 
    1/19 Sin[1 - (1866 t)/11] - 1/22 Sin[3/13 - (1885 t)/12] - 
    1/12 Sin[3/7 - (4191 t)/29] - 1/25 Sin[11/12 - (377 t)/3] - 
    1/7 Sin[11/9 - (1552 t)/19] - 1/4 Sin[5/9 - (377 t)/5] - 
    4/15 Sin[1/9 - (509 t)/9] - 2/11 Sin[3/10 - (377 t)/12] - 
    9/7 Sin[7/6 - (201 t)/8] - 24/25 Sin[1/8 - (88 t)/7] + 
    82/15 Sin[51/13 + (44 t)/7] + 7/9 Sin[4 + (132 t)/7] + 
    3/11 Sin[34/9 + (377 t)/10] + 1/9 Sin[75/19 + (1627 t)/37] + 
    2/7 Sin[2/9 + (553 t)/11] + 1/8 Sin[29/7 + (377 t)/6] + 
    3/11 Sin[25/7 + (553 t)/8] + 1/14 Sin[13/3 + (2023 t)/23] + 
    1/11 Sin[27/11 + (377 t)/4] + 1/33 Sin[9/7 + (1307 t)/13] + 
    1/10 Sin[55/12 + (1175 t)/11] + 1/42 Sin[47/11 + (1131 t)/10] + 
    1/28 Sin[71/18 + (955 t)/8] + 1/13 Sin[11/5 + (2111 t)/16] + 
    1/22 Sin[40/11 + (1244 t)/9] + 1/17 Sin[35/18 + (754 t)/5] + 
    1/18 Sin[49/16 + (1797 t)/11] + 1/13 Sin[21/8 + (1640 t)/9] + 
    1/17 Sin[7/4 + (377 t)/2] + 1/37 Sin[5/8 + (1753 t)/9] + 
    1/50 Sin[42/11 + (1709 t)/8] + 1/28 Sin[23/5 + (2419 t)/11] + 
    1/60 Sin[5/4 + (3104 t)/13] + 1/70 Sin[49/15 + (5391 t)/22] + 
    1/36 Sin[25/6 + (754 t)/3] + 1/24 Sin[8/5 + (3349 t)/13] + 
    1/99 Sin[7/3 + (2972 t)/11] + 1/63 Sin[159/79 + (1885 t)/6] + 
    1/64 Sin[16/15 + (2375 t)/7] + 
    1/75 Sin[25/11 + (2463 t)/7]}, {13.1114 - 
    1/97 Sin[10/9 - (3041 t)/11] - 1/46 Sin[9/8 - (2375 t)/9] - 
    1/28 Sin[1/20 - (3952 t)/17] - 1/27 Sin[2/5 - (3217 t)/16] - 
    1/51 Sin[8/7 - (1753 t)/9] - 1/21 Sin[11/9 - (2287 t)/13] - 
    1/19 Sin[3/7 - (1866 t)/11] - 1/21 Sin[3/5 - (1885 t)/12] - 
    1/15 Sin[7/10 - (4191 t)/29] - 1/16 Sin[4/5 - (377 t)/3] - 
    1/12 Sin[16/13 - (1175 t)/11] - 1/34 Sin[4/13 - (1307 t)/13] - 
    1/11 Sin[4/13 - (377 t)/6] - 1/6 Sin[1/11 - (509 t)/9] - 
    11/23 Sin[5/16 - (201 t)/8] - 57/8 Sin[22/23 - (44 t)/7] + 
    27/13 Sin[47/10 + (88 t)/7] + 3/4 Sin[2/5 + (132 t)/7] + 
    6/11 Sin[25/7 + (377 t)/12] + 1/76 Sin[14/5 + (377 t)/10] + 
    1/8 Sin[70/23 + (1627 t)/37] + 2/9 Sin[13/3 + (553 t)/11] + 
    1/6 Sin[25/7 + (553 t)/8] + 1/8 Sin[1/4 + (377 t)/5] + 
    1/25 Sin[49/12 + (1307 t)/16] + 1/20 Sin[68/15 + (2023 t)/23] + 
    1/11 Sin[8/3 + (377 t)/4] + 1/46 Sin[10/7 + (1131 t)/10] + 
    1/48 Sin[41/9 + (955 t)/8] + 1/15 Sin[23/11 + (2111 t)/16] + 
    1/20 Sin[11/3 + (1244 t)/9] + 1/15 Sin[37/16 + (754 t)/5] + 
    1/26 Sin[19/6 + (1797 t)/11] + 1/16 Sin[34/13 + (1640 t)/9] + 
    1/29 Sin[9/5 + (377 t)/2] + 1/43 Sin[17/4 + (1709 t)/8] + 
    1/35 Sin[19/5 + (2419 t)/11] + 1/36 Sin[9/5 + (1131 t)/5] + 
    1/65 Sin[7/6 + (3104 t)/13] + 1/43 Sin[27/7 + (5391 t)/22] + 
    1/27 Sin[37/8 + (754 t)/3] + 1/32 Sin[18/11 + (3349 t)/13] + 
    1/68 Sin[31/12 + (2972 t)/11] + 1/79 Sin[1/10 + (1131 t)/4] + 
    1/72 Sin[3/8 + (3594 t)/11]}}, {{287.833 - 
    1/39 Sin[8/15 - (1552 t)/19] - 1/36 Sin[9/11 - (377 t)/5] - 
    1/22 Sin[11/8 - (553 t)/8] - 1/32 Sin[5/9 - (377 t)/6] - 
    1/5 Sin[7/11 - (1891 t)/43] - 1/16 Sin[5/9 - (377 t)/12] - 
    1/10 Sin[11/10 - (201 t)/8] + 19/7 Sin[19/6 + (44 t)/7] + 
    2/5 Sin[1/21 + (88 t)/7] + 1/8 Sin[2/5 + (132 t)/7] + 
    1/9 Sin[1/72 + (377 t)/10] + 1/37 Sin[27/7 + (553 t)/11] + 
    1/9 Sin[24/11 + (509 t)/9] + 1/31 Sin[16/7 + (2023 t)/23] + 
    1/60 Sin[11/5 + (377 t)/4] + 
    1/43 Sin[9/5 + (1175 t)/11]}, {19.1795 - 1/39 Cos[(1307 t)/16] - 
    1/21 Sin[9/13 - (553 t)/8] - 1/37 Sin[1/11 - (377 t)/6] - 
    1/11 Sin[14/9 - (1891 t)/43] - 1/21 Sin[8/15 - (377 t)/10] + 
    23/9 Sin[51/11 + (44 t)/7] + 1/3 Sin[1/5 + (88 t)/7] + 
    2/7 Sin[2/5 + (132 t)/7] + 1/12 Sin[79/17 + (201 t)/8] + 
    1/7 Sin[3/4 + (377 t)/12] + 1/15 Sin[9/2 + (553 t)/11] + 
    1/14 Sin[31/9 + (509 t)/9] + 1/27 Sin[55/12 + (377 t)/5] + 
    1/39 Sin[13/5 + (2023 t)/23] + 1/47 Sin[5/2 + (1175 t)/11] + 
    1/72 Sin[51/11 + (1131 t)/10]}}}

It goes without saying but these equations are very long. And on top of that for the integrand you have to square these equations and multiply by the derivative. That's a pretty complicated integral, which on top of the the large number of terms, I imagine, is also highly oscillating. Since the integrals are so complicated, I imagine that's what is causing the centroids to be inaccurate at times. Whenever I run my code, I've also been getting an error message, like this, which seems to support this.

NIntegrate::levmaxord: (287.833 -1/39 Sin[8/15-(1552 t)/19]-1/36 Sin[9/11-(377 t)/5]-1/22 Sin[11/8-(553 t)/8]-1/32 Sin[5/9-(377 t)/6]-1/5 Sin[7/11-(1891 t)/43]-1/16 Sin[5/9-(377 t)/12]-1/10 Sin[11/10-(201 t)/8]+19/7 Sin[19/6+(44 t)/7]+2/5 Sin[1/21+(88 t)/7]+1/8 Sin[2/5+(132 t)/7]+1/9 Sin[1/72+(377 t)/10]+1/37 Sin[27/7+(553 t)/11]+1/9 Sin[24/11+(509 t)/9]+1/31 Sin[16/7+(2023 t)/23]+1/60 Sin[11/5+(377 t)/4]+1/43 Sin[9/5+(1175 t)/11])^2 is a Levin function of differential order 1088 which exceeds value of option "MaxOrder" -> 50. Treating (287.833 -1/39 Sin[8/15-(1552 t)/19]-1/36 Sin[9/11-(377 t)/5]-1/22 Sin[11/8-(553 t)/8]-1/32 Sin[5/9-(377 t)/6]-1/5 Sin[7/11-(1891 t)/43]-1/16 Sin[5/9-(377 t)/12]-1/10 Sin[11/10-(201 t)/8]+19/7 Sin[19/6+(44 t)/7]+2/5 Sin[1/21+(88 t)/7]+1/8 Sin[2/5+(132 t)/7]+1/9 Sin[1/72+(377 t)/10]+1/37 Sin[27/7+(553 t)/11]+1/9 Sin[24/11+(509 t)/9]+1/31 Sin[16/7+(2023 t)/23]+1/60 Sin[11/5+(377 t)/4]+1/43 Sin[9/5+(1175 t)/11])^2 as a non-Levin function.

I've tried a bunch of different integration methods, but nothing seems to work better than the method Mathematica automatically uses. I wanted to have an alternative way to calculate the centroid, so I could compare my results and see how far off some of the more reasonable answers are. So I came up with an alternative method that turns the parametric plot into a region using DiscretizeGraphics and then calculates the centroid using RegionCentroid and I believe that this gives more accurate results. All three x coordinates for each method are different by at least a couple of pixels. This may not seem significant, but based on the size of the image 1 pixel in the x-direction is approximately 1 degree in longitude, which could make a significant difference when I'm calculating the slope. Do you know any way for these integrals to converge more accurately or perhaps some other unforeseen problems that are causing the numbers to be off? Any help would be greatly appreciated. Sorry in advance if my question is obvious. I would still consider myself a novice when it comes to Mathematica, so there's a lot I don't know still. Also, sorry that my post was so long, I may have wrote more than I needed but I wanted to give enough background to make clear the full context of the problem.

$\endgroup$
3
  • $\begingroup$ My first guess (I don’t have a computer at hand) is that for the third figure, $x^2 \approx 84000$ and rationalizing the coefficient 2. Pi n to with 0.01 introduces small errors in the integral that when multiplied by $x^2$ are significant. Why not use the exact 2 Pi n for the coefficient? $\endgroup$ – Michael E2 Jun 14 '20 at 14:28
  • $\begingroup$ The method here could probably be adapted to do your integrals exactly and efficiently. You’d need to add rules for the integrals of cubic products of sine and cosine. $\endgroup$ – Michael E2 Jun 14 '20 at 14:34
  • $\begingroup$ Actually I don’t see the need to Rationalize[] at all. It’s cute for the “person curves.” But it adds an error, which again is fine for person curves. $\endgroup$ – Michael E2 Jun 14 '20 at 14:40
3
$\begingroup$

There's an easy morphological way to get the centroids:

img = Import["https://i.stack.imgur.com/VrChT.jpg"];
ifcs = ImageForestingComponents[img];
centroids = {2, 3, 4} /. ComponentMeasurements[ifcs, "Centroid"]; 
(* note component 1 is the background, we only care about components 2,3,4 *)

Show[Colorize[ifcs], Graphics[{PointSize[Medium], Point[centroids]}]]

coronal centroids

$\endgroup$
2
  • $\begingroup$ Thanks, for the info but this doesn't really answer my question. It's a good alternative method, but it doesn't really answer my question. I was wondering if there is a way to make the integrals give more accurate results. I've spent a long time on this method in particular and it would be unfortunate if I have to abandon it entirely. I realize now there are easier ways to do it on Mathemtatica but I'm thinking that this method might be used in other programs at some point, unless there are similar build in routines, like the ones used here. $\endgroup$ – Jacob Jun 14 '20 at 19:30
  • $\begingroup$ I've spent a long time on this method in particular and it would be unfortunate if I have to abandon it entirely - this is known as a sunken cost - it is not unfortunate to abandon something you worked hard at if a better solution exists - in fact it is quite brave. I hate to be the one to say it but your method is really painfully complicated even though it shows you are clearly very skilled with the parametric curves formulation. You also posted an extremely long question and I took a lot of time to decipher that your primary goal is to get the centroids, which is what I've done. $\endgroup$ – flinty Jun 14 '20 at 19:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.