circle problem, circumscribed to another circle and the area, area per image detection

5 views (last 30 days)
I have been working on a transportation problem, I have a series of coordinates that always vary, but for this question let's suppose that I have 13 coordinates, each chord represents the center of a circle with a radius of 400m.
the circles are on an x, y plane. The idea is to find the total area that is inside each circle, the problem is that it is not simply multiplying the number of circles by (PI * R ^ 2), because many of them have 1 or more shared areas that do not they can be added more than once.
These circles may or may not intercept each other, through the program that I am designing, I even obtain the distance matrix between each one, at first I thought that the best way to solve the problem was mathematically, but today I realized that it is a problem of mathematical complexity .
So looking on the internet I found that more people have had a similar problem and that it can be solved analytically.
It turns out that there are people who say that the problem could be solved if each point is plotted on a 2d plane, then a 400 meter Buffet is placed at each point, each circle is painted in one color and then the background in another color.
Finally, an algorithm is required that is capable of finding the area that is colored, suppose that what is inside each circle is yellow and the outside is green.
could you please help me with this
attached document excel file
1)Once I have plotted the figure, how do I put the 400 meter buffer at each point, I know the x, y coordinate of each point)
2)Supposing that I manage to obtain a graph with the points and their respective buffer of 400 meters of radius, how I paint everything that is inside the buffer of the same color
3)how do i find the colored area inside each buffer
T=readtable('mini2.xlsx')
H = table(T.lon,T.lat);
plot(T.lon,T.lat,'ko*400')
plot It has no buffer, I don't know how to implement it
Thanks for your time.
examples:
  3 Comments
Andrew
Andrew on 26 Aug 2021
Edited: Andrew on 26 Aug 2021
Are approximate answers OK? Since the radii are always 400 m, I would consider finding out ahead of time how much two circles overlap as a function of distance. At x/r~=0, the second circle has basically zero area, and at x/r~=2, it's basically a whole circle. Then you can simply pairwise see which ones overlap and apply the function. The only serious issue I see is that if you have three or more overlapping, you'll overcount the overlap, so I suppose this would work best if it's a pretty sparse map or the little bit in the midst of three circles doesn't matter much.

Sign in to comment.

Accepted Answer

DGM
DGM on 23 Jul 2021
Edited: DGM on 23 Jul 2021
This could be improved or elaborated, but it's a simple start:
% let's assume that 1px = 1 square meter
r = 400;
C = rand(10,2)*5000 + 1000; % test locations (x,y)
% get padded image size
p = 10; % extra padding
s = ceil([(max(C(:,2))-min(C(:,2))),(max(C(:,1))-min(C(:,1)))] + r*2 + p*2);
C = C - [min(C(:,1)) min(C(:,2))] + r + p; % shift image
% draw all circles
m = false(s);
for c = 1:size(C,1)
thisroi = images.roi.Circle('center',C(c,:),'radius',r);
m = m | createMask(thisroi,m);
end
totalarea = sum(m(:)) % total number of white pixels = area
totalarea = 4914749
imshow(m)
There is also bwarea(), but you might want to read the docs to make sure it's what you want. It's not the same as a simple sum.
Of course, if there are many such circles and they are very spread out, the image may be very large. It might be rather slow to do it by drawing circles like this. This is a different approach for the same thing, which should be faster than explicit circle drawing.
% same stuff
r = 400;
C = rand(10,2)*5000 + 1000; % test locations (x,y)
% same stuff
p = 10;
s = ceil([(max(C(:,2))-min(C(:,2))),(max(C(:,1))-min(C(:,1)))] + r*2 + p*2);
C = round(C - [min(C(:,1)) min(C(:,2))] + r + p);
% use distance xform
m = false(s);
m(sub2ind(s,C(:,2),C(:,1))) = 1;
m = bwdist(m)<=r;
totalarea = sum(m(:)) % total number of white pixels = area
totalarea = 4164733
imshow(m)
  12 Comments
David Alejandro Ramirez Cajigas
I have posted the question right here at the end, it will be very kind if you could help me
Thank you
down in this forum, is the question
Thank you

Sign in to comment.

More Answers (2)

David Alejandro Ramirez Cajigas
for Java I found this on the web could matlab be adapted?
author : http://jsfiddle.net/broofa/WRpmW/5/
// Get HTML canvas element (and context) to draw into
var canvas = document.getElementById('canvas');
var ctx = canvas.getContext('2d');
// Lil' circle drawing utility
function circle(x,y,r) {
ctx.beginPath();
ctx.arc(x, y, r, 0, Math.PI*2);
ctx.fill();
}
// Clear canvas (to black)
ctx.fillStyle = 'black';
ctx.fillRect(0, 0, canvas.width, canvas.height);
// Fill shape (in white)
ctx.fillStyle = 'white';
circle(40, 50, 40);
circle(40, 10, 10);
circle(25, 15, 12);
circle(35, 90, 10);
// Get bitmap data
var id = ctx.getImageData(0, 0, canvas.width, canvas.height);
var pixels = id.data; // Flat array of RGBA bytes
// Determine area by counting the white pixels
for (var i = 0, area = 0; i < pixels.length; i += 4) {
area += pixels[i]; // Red channel (same as green and blue channels)
}
// Normalize by the max white value of 255
area /= 255;
// Output result
document.getElementById('result').innerHTML = area.toFixed(2);

David Alejandro Ramirez Cajigas
Edited: David Alejandro Ramirez Cajigas on 17 Aug 2021
Will it be possible to make the plotted image of black and white circles to measure the area, be drawn within specific margins? that is to say that it is drawn, for example, inside a straight line and where a circle touches the "wall" of the rectangle, is that area not counted?
I could have the coordinates of the rectangle and from there I would already have the "wall" where the circuit cannot advance.I made a small drawing and I will post it within the same link as the question already solved. with an exampleThank you.
1) the point cloud in longitude and latidut is
lat lon
3.4310571 -76.5249196
3.4330517 -76.5281573
3.4329578 -76.5278992
3.4259733 -76.5169571
3.4141012 -76.5322045
3.4132041 -76.5212801
3.4147463 -76.5212291
3.4200984 -76.5090573
3.4193326 -76.5304345
3.4191612 -76.5226856
3.4195618 -76.5180319
3.4194799 -76.5200061
3.4197733 -76.514478
3.4199018 -76.5121552
3.4191972 -76.5076382
3.4186178 -76.5074345
3.4172711 -76.5061631
3.4174853 -76.5059834
3.4131446 -76.5021961
3.413158 -76.5019253
3.410971 -76.5001202
3.4120227 -76.4916586
3.4146176 -76.4889863
3.4187462 -76.5249977
3.4121893 -76.5258832
3.4266255 -76.5285806
3.4291025 -76.5167687
3.4347968 -76.5141764
3.4321028 -76.5101736
3.4271533 -76.5053988
3.4185984 -76.4869869
3.4218167 -76.4909459
3.4248689 -76.4945159
3.4282086 -76.5327531
3.4354939 -76.5187478
3.4334187 -76.5101116
3.4352097 -76.5118256
3.4277076 -76.5368692
3.4309797 -76.5358927
3.4333837 -76.5351849
3.4349393 -76.5347021
3.4327974 -76.5364751
3.4182723 -76.5367191
3.4090406 -76.5355765
3.4077261 -76.5315237
3.4109503 -76.4893841
3.4086604 -76.5217632
3.4107276 -76.5215164
3.4166607 -76.5211731
3.4197919 -76.520905
3.4226876 -76.5205857
3.4247333 -76.5203818
3.4252902 -76.5180966
3.4295177 -76.5165489
3.4295177 -76.5165489
3.431641 -76.5163424
3.434254 -76.5147331
3.4347146 -76.5147384
3.4324815 -76.5163367
3.4213276 -76.517367
3.4175283 -76.5222486
3.412915 -76.5227287
3.410417 -76.523005
3.4071183 -76.5233778
3.4253297 -76.5015155
3.4339325 -76.4964618
3.4309045 -76.4980119
3.4244061 -76.5032478
3.4351374 -76.5232265
3.4346636 -76.5230613
3.4352634 -76.5206634
3.4330464 -76.5264623
3.4329822 -76.525545
3.432503 -76.5272398
3.431776 -76.528092
3.4305833 -76.5299112
3.4284118 -76.5320919
3.42642 -76.5337575
3.4241897 -76.533956
3.4235712 -76.5346507
3.4203582 -76.5356673
3.4178333 -76.5365309
3.4186983 -76.5357638
3.4215123 -76.5348599
3.4260317 -76.5333472
3.431708 -76.53415
3.4312579 -76.5281195
3.4321469 -76.5268189
3.4104115 -76.5293914
3.4129524 -76.5282836
3.4116382 -76.5231909
3.4225161 -76.5052968
3.4223848 -76.5046797
3.4201063 -76.5077996
3.4198948 -76.5072551
3.4282383 -76.4985591
3.4307088 -76.4971941
3.4295397 -76.498683
3.4341627 -76.4954614
3.4186526 -76.5283958
3.420709 -76.5345703
3.4148539 -76.5343228
3.4143741 -76.5182786
3.4143097 -76.5204408
3.4127084 -76.5180348
3.4121921 -76.5152421
3.4101545 -76.5174233
3.4144032 -76.5129681
3.4171883 -76.5100628
3.4339407 -76.53169
3.4313357 -76.5324383
3.435523 -76.5321513
3.4296035 -76.527203
3.4293248 -76.5276801
3.4261009 -76.5289884
3.4230866 -76.5299733
3.4235895 -76.5295651
3.4198576 -76.5310217
3.4198307 -76.5307561
3.4168025 -76.5317512
3.416891 -76.5319739
3.4140367 -76.5326686
3.4140957 -76.5328912
3.4119804 -76.5336127
3.4118654 -76.5333713
3.4093833 -76.5341652
3.4096271 -76.534353
3.4072306 -76.5351255
3.4085317 -76.5191425
3.4089332 -76.5194188
3.4076618 -76.5207116
3.4112626 -76.5170021
3.4133751 -76.5148108
3.4176403 -76.510361
3.4163313 -76.5116723
3.4153969 -76.5043551
3.4084892 -76.4974891
3.4101436 -76.4990257
3.4144384 -76.503151
3.4351054 -76.5284981
3.4261281 -76.5275993
3.4261523 -76.527849
3.4229822 -76.5274654
3.4228698 -76.5277042
3.4193409 -76.5272482
3.4195142 -76.5274917
3.4174159 -76.5276934
3.4175039 -76.5274625
3.4152952 -76.5282755
3.4153355 -76.5280341
3.4136218 -76.5285115
3.4128266 -76.5289863
3.4105909 -76.5295951
3.4122988 -76.4918565
3.4172202 -76.4868609
3.414781 -76.4891394
3.4114073 -76.4938627
3.4161733 -76.4946323
3.413493 -76.4958315
3.4153673 -76.4960917
3.409121 -76.5029153
3.4107622 -76.5010752
3.4078625 -76.5035831
3.4223742 -76.5071477
3.4225534 -76.5073462
3.4296915 -76.5034007
3.4300717 -76.5034999
3.427595 -76.504774
3.4321119 -76.5024404
3.4273729 -76.5045942
3.4354533 -76.5004663
3.4320209 -76.5022017
3.4252149 -76.5060265
3.4250836 -76.5057395
3.4355899 -76.5261781
3.4306315 -76.5255236
3.4283743 -76.5244554
3.4276114 -76.5261325
3.4235978 -76.5259367
3.4216031 -76.5252152
3.4217371 -76.5227609
3.4199565 -76.5219831
3.4077286 -76.5265107
3.4108104 -76.525537
3.42046 -76.5240967
3.423772 -76.5242335
3.4260048 -76.5243488
3.4275149 -76.5116888
3.4254534 -76.5118604
3.42231 -76.5121072
3.4119459 -76.5144917
3.4120689 -76.5143116
3.410232 -76.5128663
3.4103927 -76.5127
3.4080554 -76.5104764
3.4079321 -76.5106856
3.410584 -76.5180612
3.4136912 -76.5177206
3.4146069 -76.519405
3.4281973 -76.5199898
3.4324517 -76.5213684
3.4319432 -76.5237341
3.4256217 -76.5249679
3.4257289 -76.5229053
3.4234424 -76.5215964
3.4102992 -76.5182481
3.4148131 -76.5178014
3.4211585 -76.517155
3.414896 -76.5176163
3.4173218 -76.5175815
3.4173192 -76.517375
3.4198733 -76.5172382
3.4197957 -76.5174501
3.4226928 -76.517037
3.4080927 -76.5165086
3.4072227 -76.5159668
3.4076885 -76.5163987
3.4144512 -76.4930341
3.4129628 -76.4915964
3.4176748 -76.4959094
3.419292 -76.4974544
3.4224996 -76.5005174
3.4207702 -76.4988491
3.4227033 -76.5040283
3.4205988 -76.5019228
3.4180153 -76.4994632
3.4162721 -76.4980081
3.4352152 -76.5141913
3.4335874 -76.5092426
3.4330037 -76.5074616
3.433582 -76.5052032
3.4336247 -76.5059327
3.4339756 -76.5031808
3.4327895 -76.501204
3.4337159 -76.5030869
3.4328028 -76.5008097
3.4318442 -76.4993533
3.4314239 -76.4980551
3.4344601 -76.5118336
3.4192084 -76.5207168
3.4193369 -76.5181929
3.4194627 -76.5159237
3.4196073 -76.5133407
3.4197599 -76.5103903
3.428547 -76.5330845
3.4121579 -76.4873705
3.4297281 -76.5109285
3.4261112 -76.5168337
3.4229298 -76.5171669
2) the result image is
3) the idea is that the points are plotted in a "rectangular" range demarcated by coordinates
for example:
lat lon
3.4356 -76.5369 %lan and lot 1
3.4071 -76.4868 %lan and lot 2
How do I get it to plot and measure me within the geometric limits of the rectangle, formed by the given coordinates?
Thank you
  3 Comments
DGM
DGM on 19 Aug 2021
Oof. Sorry about that. Sometimes site notifications get buried, and I almost never check email anymore. Anyway...
% This is only part of the dataset!
lon = [-76.5249196 -76.5281573 -76.5278992 -76.5169571 -76.5322045 -76.5212801 -76.5212291 -76.5090573 -76.5304345 -76.5226856 -76.5180319 -76.5200061 -76.514478 -76.5121552 -76.5076382 -76.5074345 -76.5061631 -76.5059834 -76.5021961 -76.5019253 -76.5001202 -76.4916586 -76.4889863 -76.5249977 -76.5258832 -76.5285806 -76.5167687 -76.5141764 -76.5101736 -76.5053988 -76.4869869 -76.4909459 -76.4945159 -76.5327531 -76.5187478 -76.5101116 -76.5118256 -76.5368692 -76.5358927 -76.5351849 -76.5347021 -76.5364751 -76.5367191 -76.5355765 -76.5315237 -76.4893841 -76.5217632 -76.5215164 -76.5211731 -76.520905 -76.5205857 -76.5203818 -76.5180966 -76.5165489 -76.5165489 -76.5163424 -76.5147331 -76.5147384 -76.5163367 -76.517367 -76.5222486 -76.5227287 -76.523005 -76.5233778 -76.5015155 -76.4964618 -76.4980119 -76.5032478 -76.5232265 -76.5230613 -76.5206634 -76.5264623 -76.525545 -76.5272398 -76.528092 -76.5299112 -76.5320919 -76.5337575 -76.533956 -76.5346507 -76.5356673 -76.5365309 -76.5357638 -76.5348599 -76.5333472 -76.53415 -76.5281195 -76.5268189 -76.5293914 -76.5282836];
lat = [3.4310571 3.4330517 3.4329578 3.4259733 3.4141012 3.4132041 3.4147463 3.4200984 3.4193326 3.4191612 3.4195618 3.4194799 3.4197733 3.4199018 3.4191972 3.4186178 3.4172711 3.4174853 3.4131446 3.413158 3.410971 3.4120227 3.4146176 3.4187462 3.4121893 3.4266255 3.4291025 3.4347968 3.4321028 3.4271533 3.4185984 3.4218167 3.4248689 3.4282086 3.4354939 3.4334187 3.4352097 3.4277076 3.4309797 3.4333837 3.4349393 3.4327974 3.4182723 3.4090406 3.4077261 3.4109503 3.4086604 3.4107276 3.4166607 3.4197919 3.4226876 3.4247333 3.4252902 3.4295177 3.4295177 3.431641 3.434254 3.4347146 3.4324815 3.4213276 3.4175283 3.412915 3.410417 3.4071183 3.4253297 3.4339325 3.4309045 3.4244061 3.4351374 3.4346636 3.4352634 3.4330464 3.4329822 3.432503 3.431776 3.4305833 3.4284118 3.42642 3.4241897 3.4235712 3.4203582 3.4178333 3.4186983 3.4215123 3.4260317 3.431708 3.4312579 3.4321469 3.4104115 3.4129524];
lonlat = table(lon.',lat.');
% boundary
box = [-76.5369 3.4356; -76.4868 3.4071]; % lon/lat
% process table & box; convert to meters
lla = [fliplr(table2array(lonlat)) zeros(size(lonlat,1),1)];
position = lla2flat(lla,[min(lla(:,1)) min(lla(:,2))],90,0);
boxpos = lla2flat([fliplr(box) [0; 0]],[min(lla(:,1)) min(lla(:,2))],90,0);
boxpos = boxpos(:,1:2);
% same as before
r = 400;
C = position(:,1:2); % locations (x,y)
p = 10; % padding
s = ceil([(max(C(:,2))-min(C(:,2))),(max(C(:,1))-min(C(:,1)))] + r*2 + p*2);
refpoint = [min(C(:,1)) min(C(:,2))];
C = round(C - refpoint + r + p);
corners = round(boxpos - refpoint + r + p);
m = false(s);
m(sub2ind(s,C(:,2),C(:,1))) = 1;
m = bwdist(m)<=r;
m(:,[1:corners(1,1) corners(2,1):end]) = 0;
m([1:corners(1,2) corners(2,2):end],:) = 0;
size(C,1)*round(pi*r^2)
totalarea = sum(m(:)) % total number of white pixels = area
imshow(m); hold on
plot(corners([1 2 2 1 1],1),corners([1 1 2 2 1],2),':')

Sign in to comment.

Products


Release

R2021a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!