From 221bc62e744afdb9528a5babe4664dc33dd5e7db Mon Sep 17 00:00:00 2001 From: Olly Date: Thu, 21 Sep 2023 23:27:31 +0100 Subject: [PATCH] Add Circularity methods --- .../imports/simba/simba.import_circle.pas | 24 ++++++++++++++ .../script/imports/simba/simba.import_tpa.pas | 12 +++++++ Source/simba.circle.pas | 23 ++++++++++++- Source/simba.tpa.pas | 7 ++++ Tests/tpa_minareacircle.simba | 32 +++++++++++++++++++ 5 files changed, 97 insertions(+), 1 deletion(-) create mode 100644 Tests/tpa_minareacircle.simba diff --git a/Source/script/imports/simba/simba.import_circle.pas b/Source/script/imports/simba/simba.import_circle.pas index 5146ff8e5..9cc6e4212 100644 --- a/Source/script/imports/simba/simba.import_circle.pas +++ b/Source/script/imports/simba/simba.import_circle.pas @@ -123,6 +123,28 @@ procedure _LapeCircle_Circumference(const Params: PParamArray; const Result: Poi PDouble(Result)^ := PCircle(Params^[0])^.Circumference(); end; +(* +TCircle.Center +~~~~~~~~~~~~~~ +> function TCircle.Center: TPoint; + +Returns the center point of the circle. +*) +procedure _LapeCircle_Center(const Params: PParamArray; const Result: Pointer); LAPE_WRAPPER_CALLING_CONV +begin + PPoint(Result)^ := PCircle(Params^[0])^.Center(); +end; + +(* +TCircle.Circularity +~~~~~~~~~~~~~~~~~~~ +> function TCircle.Circularity(TPA: TPointArray): Double; +*) +procedure _LapeCircle_Circularity(const Params: PParamArray; const Result: Pointer); LAPE_WRAPPER_CALLING_CONV +begin + PDouble(Result)^ := PCircle(Params^[0])^.Circularity(PPointArray(Params^[1])^); +end; + (* TCircle.Area ~~~~~~~~~~~~ @@ -213,6 +235,8 @@ procedure ImportCircle(Compiler: TSimbaScript_Compiler); addGlobalFunc('function TCircle.RandomPoint: TPoint', @_LapeCircle_RandomPoint); addGlobalFunc('function TCircle.RandomPointCenter: TPoint', @_LapeCircle_RandomPointCenter); + addGlobalFunc('function TCircle.Center: TPoint', @_LapeCircle_Center); + addGlobalFunc('function TCircle.Circularity(TPA: TPointArray): Double', @_LapeCircle_Circularity); addGlobalFunc('function TCircle.Circumference: Double', @_LapeCircle_Circumference); addGlobalFunc('function TCircle.Area: Double', @_LapeCircle_Area); addGlobalFunc('function TCircle.Expand(Amount: Integer): TCircle', @_LapeCircle_Expand); diff --git a/Source/script/imports/simba/simba.import_tpa.pas b/Source/script/imports/simba/simba.import_tpa.pas index 6e44b37d4..0af12ca3d 100644 --- a/Source/script/imports/simba/simba.import_tpa.pas +++ b/Source/script/imports/simba/simba.import_tpa.pas @@ -744,6 +744,16 @@ procedure _LapeTPADistanceTransform(const Params: PParamArray; const Result: Poi PSingleMatrix(Result)^ := PPointArray(Params^[0])^.DistanceTransform(); end; +(* +TPointArray.Circularity +~~~~~~~~~~~~~~~~~~~~~~~ +> function TPointArray.Circularity: Double; +*) +procedure _LapeTPACircularity(const Params: PParamArray; const Result: Pointer); LAPE_WRAPPER_CALLING_CONV +begin + PDouble(Result)^ := PPointArray(Params^[0])^.Circularity(); +end; + procedure ImportTPA(Compiler: TSimbaScript_Compiler); begin with Compiler do @@ -837,6 +847,8 @@ procedure ImportTPA(Compiler: TSimbaScript_Compiler); addGlobalFunc('function TPointArray.DistanceTransform: TSingleMatrix;', @_LapeTPADistanceTransform); + addGlobalFunc('function TPointArray.Circularity: Double;', @_LapeTPACircularity); + ImportingSection := ''; end; end; diff --git a/Source/simba.circle.pas b/Source/simba.circle.pas index 20f5fa4d0..fa471f37e 100644 --- a/Source/simba.circle.pas +++ b/Source/simba.circle.pas @@ -43,6 +43,7 @@ TCircle = record function Exclude(Points: TPointArray): TPointArray; function RandomPoint: TPoint; function RandomPointCenter: TPoint; + function Circularity(TPA: TPointArray): Double; end; operator in(const P: TPoint; const Circle: TCircle): Boolean; @@ -51,7 +52,7 @@ implementation uses Math, - simba.math, simba.tpa, simba.random, simba.overallocatearray; + simba.math, simba.tpa, simba.random, simba.overallocatearray, simba.geometry; class function TCircleHelper.Create(AX, AY: Integer; ARadius: Integer): TCircle; begin @@ -167,6 +168,26 @@ function TCircleHelper.RandomPointCenter: TPoint; Result.Y := Center.X + Round(R * SinValue); end; +function TCircleHelper.Circularity(TPA: TPointArray): Double; +var + I: Integer; + Smallest, Test: Double; + Hull: TPointArray; +begin + Hull := TPA.ConvexHull(); + if Length(Hull) <= 1 then + Exit(0); + + Smallest := $FFFFFF; + for I := 0 to High(Hull) do + begin + Test := TSimbaGeometry.DistToLine(Self.Center(), Hull[I], Hull[(I+1) mod Length(Hull)]); + if (Test < Smallest) then + Smallest := Test; + end; + Result := Sqr(Smallest) / Sqr(Self.Radius); +end; + operator in(const P: TPoint; const Circle: TCircle): Boolean; begin Result := Circle.Contains(P); diff --git a/Source/simba.tpa.pas b/Source/simba.tpa.pas index 89db9f806..915fdde3e 100644 --- a/Source/simba.tpa.pas +++ b/Source/simba.tpa.pas @@ -134,6 +134,8 @@ interface function Intersection(Other: TPointArray): TPointArray; function DistanceTransform: TSingleMatrix; + + function Circularity: Double; end; implementation @@ -2269,6 +2271,11 @@ function TPointArrayHelper.DistanceTransform: TSingleMatrix; Result := Transform(data,w,h); end; +function TPointArrayHelper.Circularity: Double; +begin + Result := MinAreaCircle().Circularity(Self); +end; + (* Approximate smallest bounding circle algorithm by Slacky This formula makes 3 assumptions of outlier points that defines the circle. diff --git a/Tests/tpa_minareacircle.simba b/Tests/tpa_minareacircle.simba new file mode 100644 index 000000000..45e199936 --- /dev/null +++ b/Tests/tpa_minareacircle.simba @@ -0,0 +1,32 @@ +{$assertions on} + +const + expected: TCircleArray = [[261, 87, 84], [451, 55, 83], [69, 87, 62], [444, 178, 90], [102, 292, 105], [241, 240, 55], [348, 373, 113], [475, 311, 42], [178, 418, 34], [86, 419, 30]]; + expectedCircularity: TIntegerArray = [54,18,93,13,39,48,41,65,58,59]; + +var + img: TImage; + TPA: TPointArray; + ATPA: T2DPointArray; + I: Integer; + C: TCircle; +begin + Img := TImage.CreateFromString('IMG:eJy9WtdyJDtynZW+Qz+BDfySIvSk0B+WaVftfdM1vfd2ZsghNvacA6C6usm9Ib1obvCSyDxIAIl0AOp//vu//vPHj3/7jx8//jb+9x//+Pvf/+7+b/+e/pEkSWo61kTKtq2Zpk3L9tia3LZMCnKyQvVtZmwblIaxSUlt4u/MDoztlKQ6RKW2b2qQsaIl6Nkz1oCZ2KLCMHZsbG4q4Mw2jV1WIGfx70/Lse1lbP9GK7U3sfnLonVbdrTXppz+G3rerZovaD5xLqSMLOc15f/E/2N/cUbO1UC5CcQH6zID1kdon7OzgyJ+WtcU5RAKgr6sM+7RuqFouxZqcFDtnXVbxl2D3ZWAqTV166hgDHRo3YFxJ3YbjLl1Lal93zqse88egphAdMpd4ZhtgxY2zWHdc3vOSeXgYndcBz3rRGJN7sm4IbXxrDlsGWxc4lfcse7NuK591OyPtZ89696Na1kwHriIK4oBNrcOC6rZ38bdU/dgXUY6pp7YT+NurZOqMQx08SFNgXZhsHcOk/rFlVKDp75fz7xyle6K/RwII/PEdZFwztl8gjYx6NA1UDRosAnnZgZaHBh3wSlgFLcw19aNjTtTt9T8Bm3bXFr0dcfS1blwuwZC5sYtQctIfAdxaU61J0vrbtXzwBxbt2PcEWeOKbsjg43ZNdhEd7ci7FHrUC40504MdmmfhFSauTfPokImtm1q3StXg3VgK3epaTex7qd5FGFH+pgQlEqz1+vkN5oB1IjZbEkHo0C7E21BLWCHIe+hQsA+/jZykwQ/F+LMpeceB8EmZ5bqKMkFyUS7DxgYlQLWjEoBC5ISGgA84ajC6OLvLBLPSLynVVHjH2EZnkaLaJkrUaCOR7//hSgnxIiSahAYcKAD+Qy7MImMGiZ1IeXPSP5jFZhASDnd8xUHNkI/wAoh8E4bOqNyHX0D2/5LDm7NvTZ/XrJyauxFZrIIxJz4B8akZ5nUglqG/Ey7/Bj8CFp7kRVuiV2jn9FPsdaUfyc0c1jttthNC6HeTA45w7Y1bzLzbY7PQAhlHWnXrJbmHWMncnMFXzelbI7vGvKkvchvKvxoI2eKFBiRzrcUYKA4Te7cwGVpcK4t3R4EPvoE9p6BgqXkvgwMfn4k0NCaBuY/UbRJEDETaUqL6stmHX3ow7gdO6R6WsDCYDuI/G7BnTkVdxfoHhAjbVhOOCz7XLxtpQilnsz0AZxaKeOIer8UZCFKol1NzYDCGR9uzI3YRxgYezFEpHVU0hl3HPRPEQrG2icRnJJWj1HgNRDaSgcI3u+ekCgtLRTiPSW1jPZ7Pg/kVdRxzBbtKvVGRjfzyWVS5Tiu6a+b7xrQN4+1BWgWkftt86/AidYn7k9x/6qZhb4b3BO5CLh/3ezEJvzS5YF7Gpsb3E7kPirpbjTbId8ynsDQPffMN2vaUzXvxO1G8G2leaFAB1v4tnmpuAlv+baplIQgFprXis6wol7k+qbn3ijGwwB6sXmoZl/NW+UKBM9ebC5lf77vvVILyqxOrDeQi+CXvvmoXAUP9Np4Vj6DqfrmC+MCC7dObCIrwtp8E2Gmr5LDD/RTadXNQvOXMi9SmG/+VmaGyfvmh1I3vLIbSiPmdvirX8Insz/c2bdgI7Sy/djMWEyg6VePfUZFAlMMzYYqtZPYLFjJYOMGvtmxR4Y7NfLNnmqjG4QTLoulzJGVKmM9dQ8Xo7b6CJhNhFR1u7MLgSZIBB0zDyXurQohHzSU34zdQTgamQnrKm1bQwUsV9tBHMnDVjcRb1DW/WIR0re1YHgFgtGUVRZ03g3UQ8tSuYHhLYs+7FWBuoycXRXBBQtezA02vo3A7VlbmONIQa2h2qtQ0Dsygb1g4Y1MlEJuDZG1phLxWrwZ5raDrswBbeJVHUpF0EwTSSJDjxy8mgLxDcK0Cao4AI/z6XKdUzqO51EL+KurMAzNjelxh770v8avdxLlpAfleWBCj15Wm68YOzbHKhCeMdO0pLAkecF+ZYEy8hnsFTk1kobq9oYFJMgkG8Cf0AZzaV6CRf6F6jr3SqqSWfzkWu0a2dUr0EEoR+qQXK/gUtHq0ErEfU/LRMvL80zfF52sSUOb6cnJQyO/plS/4jcqbZZtzOTDSruF9qBsp2r3y/6mrcpgJY7tbmU41P+uXWl3K22K66lGi+2sbOferRoyiRWkA30OGAhrdieQuCGoHE6g6/mqI8ujEfMv65SuHZfgzIxJnkTvBylFC6SFIlw9WAzLlinJ28rCjXXyDiuHPW/iqEKzwPMsOoLfdjj8zFPg/5dh31fE29IdQ9Nv7lxHAa7kjGPvK6ae2hNjXxSZ1mipffbi17B9gyBT40GTweyKDr4sWQlCSxPBIReTKzrVGTwUrfTMEsnwoygSoYAdQlwRsL0Sm0LiMsAgnsAFgMMo1hehAvdgpQR3/FmeYvuKjl12KxhT+V++6tExTL9LRmBaF5Q/9SdEw8Nt2xyqRI7owjDh7oOeagDZHSfUJp7lxYldA3cA7njwJ8B1gRkkUX2tiyxUHAOlA/ghg2udsStKa8FBWAi3V7gTVeoqtb8BFkZVyB8BVYxj6DVcAUHE1bmLnN6przp5lg+QNlhLzbHlz+gn/lDA5VZkNeGOBLawPwF4DaNcQ3TB8QiM9gHEJYhHEXECjkc0Ma8WERmq6bZOrqWENk6g5EPCOyRw+yoSCsyiC+22MUaHOY7T7VYRXeimjyTakdJRRXxB9MEZwI16WCbCzxsQx/xjd4M/gKWB8OoPnoRUhAyg0jWQtukQf6+BGIMwaeJyPx0FlTXcOeb6qFGjE33yxNqPM4LXFsZHyLV+/Mf7sAMIXyNy0AP4MIxsaXlmCHTmmgNNyiiOVLqQtaA7DBFFMYXB135D2lfXz7BZ7VtTnodTDaPzpdXeZPetDxEtYA43++744mZEzqpbDt/gIhAsJtUepGN2LIQONun7iNU9TR6jPVVYDbA28Ly9I36gRa3TMeZYR7s2IF/6ZFBTk5NqfeWlUE1js1NLtUbL6hA5o5a/sDkkD8OLSuRa5x5LS1VOAfdC5P5KC/8lSDJrTK0GO7C7Sd5lVTr7QtrApbpU3UamrDKeqD/GlPx7liocclYDPHL7qAXeD2xtdHqwcF7MnDXp9iZd6kvRoyZtjCr8e1p0HeIWvDmqB0+qRfYd5jj2KaUGSJlZsoqIW2B4VKEM1nUdyxudyDM8tnS56wtOODIvdQJwe+rbDWO3gSlkiG3mbsAuLA9gKHg5jS4dqMZbii1af8ScMSNI2laANcRMAetqxBWSxRhPaaghEthil2YzWhd3qhpcoBn+HCOALYAYWOW/rxAWJzylzLyYQtdcfWWMJjGPLGOWqnj6AO1AIhfBercEtngh/sAKZld5lpF4m/u4hpK4B3+jkpWoGZeYRlSY2AMLn220ZhqUZ4AOvcTqKWBd4j1LIm7cTBIhjSVHE3upE8flapL7vOkO4D5X0Q+JsUZwQ5mrnMVBAC88mHH/REOmCLubyFtWUvMVcooQc7Qhss7LcgFnHtgQcEF/rWMKYYt5xDoKwCn697gMAneo0DqmdOGBV8aniImOR1heCaqpKjy2FwExFoIaoPAWdrICO+EDxu0KZgqV+A3s8+4K5S9crpl9h5hrR+anS7kAy1XzrsFywXrhZQBb8z1wIHkYuQ+RgoTLHd0PXTFH8zzb1gFrGMf7gulBM4X2A7s/or0n/xJZi8gMyCAwXCEJdslSvyOB1AXmPKZ/1eweDaVE6pEgIFtR4NiHOFWZgif+diIpBbeh0abyCQRP+Eup+QtY9ws8mhPc8ODRv4C1VjJrgGXUp5D+umMd2PDyapK3lF+tUOe8xmhgkxoSl+hMMMQsMjC0srE/CK2B615kXSJzFWO+oC/sMIL1hJQHMCXXJTnDVDvcr0aQ29IbUA923lQaqWkTCFXRD1GLldCC9x3nFaFEHnMjmgH2uQarSWAaYXS8hq8rqsgzTCIHqyaP59C83mj4iFIR2+Y8q2Kb6MT7xT5fHVNFi4KxuRTf4a0NV+bFp2FzC/RYcPc6VBnRf4ju8W730tZ1q5LrLJ2oiu7DaHd0RhUOB946wmJd8cHj5DAsPfqIe0Ou1Ett6dacZXyu2dYUZRP1InRstdO8XCv0qmNupCPKzYSs6RC+zwCQaQofPnr8MbeSGWfQVtE7s3zpMCuRE/tu7nQd1DCZrn54/9y1PE7znvCdoJl9M/diRtAWVLMPR51EyBxn2wfpt66LkxYKhJ4glPKbkIV9Mo+aS2MNsqX7uDuZQcFEc4e6kZI8rImt8JJmqsNzKK/F3HZjnqX/phTbCrAt3SJcYxo1SttFZH6h+leohJt64OshdxWBS0TeV12AxNktQoXD/x0YXQ4sqj32cXp/s7pEy1V+L6C3/XBNiP1YE3+ANPLT6latFsCKDvu046Hlhe1lVbCjyHpAZlbJFWbU/xaZy9KiotrcP57jLzZgmQynCRGcyL6P81slsMVZnvqHrmYAdqV4Xm6cVxV1r1SSCtYKsAGz3vYm8lGJwrtiETS0Y3k9claFPWkn03DzSIFLWnRdQ5fIPSTRZwWBRMhOQM4Z3HY3kS/+JqCoAHeFOq2i3uRtHD0LMK54KVgDnGBDb3JIupqXtCPIlN8XbIH9S1dlmRTvB+IF/cTq3f1KF6CFekc2e494rT2H06bwo3c/20RhdQtIjsGqcGjvWXuMlff+yEsIqQoa8ruIOr+OwJC8G0JI1FGBVjsAtCOtbytY8BKu45/tciVnx7u0guE+CTcdg1L6lDWt7r06uo3JITT2VJgZQPaANqPhhmVHnlTrWmMauk1jl55/h8t075f5oyoylNdqX60sfPjA6Zzo+c6xFGdxLzseMjR3wlAF/MM/PFwbzgNjdnSL70jnlY83uoYeKh3s7dHoOUQXhv7M2gz7nst5cxVsD0S1MECLZ5i6YdmtomxfDwe0Nglq6jbF0bd2FeWzIMyv4oqooZ4hGhHFyOLlJLo3FKKOkrbOmqDOCmgS2LVwSmJo73FOrPjOOcWmPyUn5VGTd9nOf6sxUmEeHnZ3IlUvmNVmX0/WQTDvpJ2ejgJ3qBfrsonso5ek0Bzrwbrsyxfd6Yo7C9fSvW+b83BJ3fd95+GOuu+5i3AnPVg1OzKltebY941NXpk8RwJLQl7lfFQJKd+nbJxNYZR+E7ho4p+ORp76E8bpqbOosYJvNmNQrvFDfbEuocx7fXijp2OCbtUtD/7By4tjf3BHjjwL1OoPEfHjpRccSCKVL+u5spAfTUfrXe+zV0aPS5F+BQGp/IRza/nX8FlgniKQJMFn/eTpXovAPYQQXpkiCvEJPkrcYwXAPtjs1ManiG35r+NgUfrMzv0XDDa+N4wZh9wdqJHCCzPdz+0Gne6pWi5fUDqKGwMx/BJu8Ds8DTzpkrYwCgHsfOK1s2eCSgYc/w/+d+optxiX7zHho69rfaihr6Ick2ZXTN4TB9K5rnn54JHx3OhV1kHjLqjrFtjAOFb3ia7PaDBNXa9Pw++k/ArtEDKhoixegfel9Tw2O+qcr7hp7MhAleq00MLPUtEyRRDwvylm3wOXIjBrdJX2+EBJnfm2fASR7xTl2KGeChJYQSpd+d+Jx3zClG4FqStWpyqNWJ7U8XshY2zIcEbaHlnxLSrVkS7l0vBopg0q1LMls5nL5Ar93QwfmfVNeD7lNvC2IUd2vkctPdbyEmX9XLJ8J6bsXHddiSrXXNk1Ckx0Cn5VnV5gnQ+o4acmCOqqYyuA+SKVVQRV28yUxhtGGz93qE/PpRBeCGTenVpGp34/F42QIX1lfAduKnifQaM87d8g+pyht7JTZhQuPDRZ+3vb+KeoMwzfgC6v0e80jDoUYrTRMw30HZ2qeEcwxiaNEd9fvumZ6mGYjsgHtLTS7xg6yvBzqX5zE/QTuwqfmTD8etcj7eocan7FOW1e6nYY1dRYjSw1lRLqKjFq/wspzbAraejZ8JIju6FJ+gLmCDo7l5BFdRXa9EZFe41gCNX2VxnHpgQUZqN/Zcd5gfG168pYJkbHu9V+9xm/93QM3uwTrXRSschp5W9+33CIbTpTj5OKg0wrzjEN9Jr8049ziB5nUTWbbllTGb7+d8vInNK1vvPgLQ1FBv93od2sBZdPbKeMHL4dd/lgTUpH3p8KTSmD8LsTQlH7m54nNM5uiGfdAO7qpidRcu76vDUK1rnqd67ouYqayxBND78gL31hleNYpmBcAyZQjjfAp0gFDYhu4qeBPF5fMwJ+BOmCT69oU6Nnex9jqvt/ppB3b3x+uZNGWlbeuQIkAKQVAHPdk7jn3K/30Fvf19wjmtOhLlTfIqzHByy3U3LoIo6fGpD4IOIlCbxW31UcftXnlNfGB/q9QLuGxjNVbr6q4eFpHzngCrwb1d1xnrdCqPNBCeB2PRifugk4kh6uIfaVH6vqE7gHowR7HMj3/OaF1V+qbHeKfboB+UFb7vTR8r6OmOeQ8IrTZUf0Wx1iPe2JE+Ha/BeO+yrBTxHdX5FzkmC4R6p6LnTa8TXYDfb9HphXGwzuMGj5QoJSre0O1AeAWCOqXrv2suivD+K9mN+MBCiCMn8ryKPXI8gfMIWU9MR/atWHwGfs+afcwJ2bUIUOEXRfWMmiiEj4FVjqL+ImyN5vVt/A0oJQnPHE1/MF7hwzeuCnnWOjhH+tIJT6zyZ2MaHfsCQUrbPwxdC1nCie/1hSDHG0vOHXulNldqnlxvhEz6Lff4yVWH1HvAUxL6F+VYF4xrXqSn0R8n+mYi1X0aFSTN+4NGFC/hPbzVJYtc29rwu7ZSkcf1RGv/jPNEahPi8LVqim0JpX1So/pi9r3liLYpZovccWC5DwJz/VYAm3+q7//+nf7Q//72//BFtXNF8='); + + TPA := Img.Finder.FindColor($C0C0C0,0); + ATPA := TPA.Cluster(15); + + for I := 0 to High(ATPA) do + begin + C := ATPA[I].MinAreaCircle(); + Assert(C = expected[i]); + Assert(Round(ATPA[I].Circularity()*100) = expectedCircularity[I]); + + //img.DrawCircle(C, 255); + //img.DrawText('%.2f'.Format([ATPA[I].Circularity()]), C.Bounds, True, 255); + end; + + //img.Show(); + + Img.Free(); +end.