Conference rusure::math

If you have a map of the object you can get a very rough idea of the
area by cutting out the piece, weighing it, and weighing a known
area of the map. The area is the ratio of the weight of the map
by the weight of the known area. The accuracy is not too good
considering, that  materials such as paper are not uniform in thickness.

The next accurate method is to use a planometer, which is a 
mechanical device that integrates the area under a curve, and there 
are models used specifically to measure the area inside a simple 
closed curve.

Finally, you can measure sets of points (X,Ytop-Ybottom) and
integrate using some method, like simpsons rule. This rule requires
an even number of points, in which the X values are all separated
by the same amount. 

If anyone knows how the planometer works let me know. thanks.

TOM

     Then there are the Monte Carlo methods.  Throw darts at the
     map in such a way that each square cm of the page has the
     same probability of the dart landing in it.  Then count the
     darts landing in the image of the plot of land on the map,
     vs the number of darts thrown ....
     
     :-) :-) :-)
     
     Dan

    By irregular, do you mean straight lines connecting various points?
    Or are you including curved perimeters as well?
    
    If the first, and the points all have rational coordinates, there
    is a fairly easy method.
    
    --  Barry

    xerox it onto graph paper with tiny squares and count the squares.
    tom

    The boundary lines are curved.
    

    Weighing works well if you use high quality paper and a good
    balance.  Counting squares also does well.
    
    Rule a large number of equally spaced parallel lines across it and
    measure the length of each inside the perimeter.  The area is 
    approximately the total length times the spacing.  This works for 
    concave and convex perimeters, but it is easier if the perimeter 
    is roughly convex.  
    
    You can do the above in several different directions to get several
    values and estimate your error.  Or you can do a sequence of
    measurements, decreasing the spacing by a factor of two each time,
    and taking the limit as spacing goes to zero.
    
    You can do something similar by choosing a center point and then
    drawing radii to the perimeter.  You end up with a lot of wedge
    shaped pieces and the area of each is approximately
    
    	area = pi * (average radius)^2 * angle / ( 2 * pi )
    
    The planometer does a mechanical integration of this formula, usually
    using a center point outside the area.  It thus adds in a bunch
    of positive values tracing the far side of the curve, and then
    subtracts while it is tracing the near side.  It uses two perpendicular
    wheels to measure radius and angle (actually a lot of chords) but
    I don't remember what it uses for the multiplication or summation.

    In an interesting follow-up, I received the following mail from Eric 
    Tatara, which I post here with his permission.
    
    --  Barry
    
    
From:	DPE::TATARA       "ERIC APO-1/C4 DTN 289-1603, (508)689-1603" 12-OCT-1989 17:12:33.16
To:	DWOVAX::YOUNG
CC:	TATARA
Subj:	Calculating the area of an area...  (MATH Note 1046.4)

Hi Barry,

I am looking for the algorithm that you mentioned in the Math Notes Conference.
My problem is that I have a bunch of points defining the corners of an area, and
I'm looking to calculate its area.  The area has no curves, and can assumed to
have no overlaps.  The angles need not be right angles.

Do you know where I could get that algorithm?

Thank you,
-- Eric

P.S. here is an example that the algorithm should be able to do:

   ------------
  /            \
 /              \
|                \
|        /\       \
|       /  \       \
-------/    \--------

    As it turns out the method misspoke somewhat in .4.  The method that I
    had in mind was really intended for use by a human using graph paper
    and did not translate nearly as well into a computer algorithim as I
    thought that it would.
    
    Somewhat abashed by having mislead Eric, I spent part of this past
    weekend tryng to come up with a simple and effcient algorithim that
    meet the criteria that he described.  I think that I have suceeded with
    a method that I think is rather elegant (ie. tricky?, unobvious?).
    
    I will post it in the next reply behind a form-feed so that those who
    want to work on it themselves can.  I found it to be a rather
    interesting puzzle.
    
    --  Barry

Eric:

	What follows is a FORTRAN routine to calculate the area enclosed
	by the line segments connecting a series of points on a 2-dimensional
	euclidian plane.  The only assumptions are that the lines never cross.

	I think that you will find this to be a fairly efficient routine.
	Also, if you are planning to convert it to Bliss, remember that
	Fortran array index's are reversed.  So that "Points(x,k)" in FORTRAN
	would be the same as "Points(k,x)" in most other languages.

!-----------------------------------------cut here-----------------------

	OPTIONS /EXTEND_SOURCE
	SUBROUTINE CALCULATE_AREA( Points, Length, Area )

    	!++
    	! This routine calculates the area enclosed by the line segments
    	! connecting a series of points in a 2-dimensional plane.  The
    	! points are passed in the POINTS(,) array.  It is assumed that
    	! the linesegments do not cross.
    	!
    	! This routine works by application of a method that I call
    	! "the Cutout Method" (Though someone else may have discovered
    	! it before me, so it may have a different official name).
    	! This basically accumulates area by calculating the area of the
    	! trapezoid underneath a line segment down to a base line (Y=0 in 
    	! this case).  As the line segments progress across the plane
    	! each adjoining trapezoid is calculated and added into the total.
    	! At some point that line segments will move back in the opposite
    	! direction to begin closing the curve.  When this happens the
    	! polarity of the math will change and the trapezoids will be
    	! calculated to have *negative* area.  These are added in just 
    	! the same which has the resultant effect of subtracting the
    	! true area of these reverse direction trapezoids.
    	!
    	! It is called the "Cutout Method" because it simulates the
    	! effects of cutting the closed curve out of a piece of paper
    	! in the plane.  As we cutout the top lines we accumulate paper
    	! below the cut in the form of adjoining trapezoids.  When we
    	! reverse direction and cutout the bottom line, we are effectively
    	! removing the paper below the line, again, in the form of
    	! adjoining trapezoids.  As it happens it does not matter how many
    	! times the lines reverse their direction (along the X axis) or
    	! how convoluted the resulting curve is.  So long as the lines
    	! never cross, the math still works out correctly (think about
    	! it).  It also does not matter what the chosen baseline is,
    	! whether it is above or below the curve or even if it runs through
    	! the curve.
    	!
    	! The one anomally is that depending on the order of the points
    	! in the passed array with repsect to the direction they progress
    	! around the curve, it is possible for the resultant area to come
    	! out negative.  This routine checks for that condition and sets
    	! it positive if this happens.
	!--

	IMPLICIT NONE
	INTEGER Length, I
	REAL Area, Points(2,length)
	INTEGER X /1/, Y /2/
	 
	If (Length .le. 1) then		! Not enough points
	    Area = 0			!  Set to 0
	    RETURN
	End if

	! Start with the area of the trapezoid connecting the 2 ends of 
    	! the list of points.
    	!
	Area = ( points( y, 1 ) + points( y, length ) ) 
     +  	* ( points( x, 1 ) - points( x, length ) ) / 2

    	! Now "cut out" all of the other trapezoids.
    	!
	Do i = 2, length

	    Area = area + ( ( points( y, i ) + points( y, i-1 ) ) 
     +			* ( points( x, 1 ) - points( x, i-1 ) ) / 2 )

	End do

    	! Check for a negative area.
    	!
	If (area .lt. 0.0)	Area = -area

	END

    Re .9:
    
    > 	    Area = area + ( ( points( y, i ) + points( y, i-1 ) ) 
    > +			* ( points( x, 1 ) - points( x, i-1 ) ) / 2 )
    
    Should that be "i"?                ^
    
    
    				-- edp

    That's a trapezoidal variation on a well-known way of calculating the area.

    Sum the areas of the triangles from the origin to each pair
    of vertices cyclically:

	area = 1/2 SUM(i = 0 to n-1) x[i]*y[i+1] - y[i]*x[i+1]

    with the subscripts interpreted modulo n.

    - Jim

    Re .11:
    
    Actually, the previous formula is better, since it has two
    addition/subtractions and one multiplication per iteration, instead of
    two multiplications and one addition/subtraction.  (Of course, the
    divide by two should be moved outside the loop.) 
    
    
    				-- edp 

    Re .10:
    
    Yes.  That SHOULD be "i" instead of "1".

    Re .12

	Neither is really optimal - I just wanted to mention another geometric
    way of looking at it.

	The following summation will be faster and also more accurate
    if the figure is not located near the origin since it avoids differencing
    the products of large nearly equal numbers (the product of differences
    is not as bad.)

	area = 1/2 * SUM (i = 1 to n-1) (x[i]-x[0])*(y[i+1]-y[i-1])

    with subscripts modulo n.

    - Jim

    Re .14:
    
    More accurate I can see.  But is
    
(14)	area = 1/2 * SUM (i = 1 to n-1) (x[i]-x[0])*(y[i+1]-y[i-1])
    
    faster than
    
(9)	area = 1/2 * SUM (i = 1 to n) (x[i]+ x[i-1])*(y[i]-y[i-1])
    
    only because of the loop shortening (n -> n-1) or is there another
    reason?
    
    --  Barry

    It does save a trip thru the loop...

    But the real reason is code generation - x[0] can stay in a register
    during the loop, saving a memory reference.

    Either loop could be reformulated to have only two fetches per
    iteration by doing two or three iterations inline and keeping
    consecutive values of x[i] and y[i] in registers.  This may not pay
    on the VAX because loads can be combined with arithmetic operations
    and this would just take an extra instruction to fetch the next value
    of x[i] or y[i], but this would win on a RISC processor like the R2000.

    - Jim

    Hokay, I can see that.  But I think that the following is closer to
    Optimal (depending on the relative cost of adds vs. multiplies, the 
    smarts of your compiler, etc.):
    
    Area = 1/2 * SUM (i = 1 to n) ( x[i] * (y[i+1]-y[i-1]) )
    
    --  Barry