Can You Calculate the Area of a Contour in Gnuplot?

Question

I've been using gnuplot for a couple of weeks now. I have large data files with 23 variables, but I select specifically x-y co-ordinate data and fluorescence intensity data for my analysis.

On of the things I would like to do is a contour plot of my fluorescing particles. I should add that this contour plot is over time so there will be several spots nearly overlapping, but this is in fact the same particle. I would like to draw contours around these spots, colour code according to intensity and have the area of the contour displayed on the graph.

I have achieved all but one of these goals for my contour plot. I cannot devise a way for gnuplot to calculate and display the area within the contour. If I could then I would have an estimate of the area of my particle. I recognise my goal may be beyond the capabilities of gnuplot, but if there were a solution then it would be very neat.

Here is my script for the contour plot which as I said gives everything I need bar the area within contours.

The co-ordinates are in nanometres and each point on the dataset is the centre of a molecule. I have taken a small range of co-ordinates because there is so much data, it would not be possible to distinguish otherwise (there are over 80 000 data points). I have also set a threshold of intensity as I only want relatively bright fluorescent particles (done with set cntrparam levels incremental 8000,5000,100000). $23 and $24 are the x and y co-ordinates respectively. $12 is the intensity.

#Contour plot of Fluorescent Particle Location with Intensity
#Gnuplot script file for plotting data in file "1002 all.txt"
reset
set dgrid3d 100,1000,1
set pm3d
set isosample 30
set xlabel 'x (nm)'
set ylabel 'y (nm)'
set contour base
set cntrparam levels incremental 8000,5000,100000
unset key
unset surface
set view map
set xrange[20000:22000]
set yrange[7000:10000]
splot "1002 all.txt" using ($23<22000 && $23>20000 ?$23 : 1/0):$24<10000 && $24>7000 ?$24 : 1/0):12 with lines
set terminal push               


set terminal png                
set output "1002_all_fluorophores_section_contour.png" # set the output filename
set terminal png size 1280,760
replot    
set output                               

Answer 1

As @Christoph says, gnuplot might not be a numerical tool, however, the calculation of a polygon area is not too complicated and can easily be done with gnuplot only. Assumption is that you have closed polygons, i.e. last point == first point, and the data of the individual polygons is separated by two empty lines.

edit: script changed to work with gnuplot 4.6.0 as well.

Data: SO28173844.dat

 1   1
 2   1
 2   2
 1   2
 1   1


 3   1
 5   4
 9   0
 8   4
 7   4
 9   8
 6   8
 4   9
 0   6
 3   1


 4   0
 5   3
 7   1
 4   0

Script: (works for gnuplot>=4.6.0, March 2012)

### calculate areas of closed polygons
reset

FILE = "SO28173844.dat"

set size ratio -1
set style fill solid 0.3
set grid x,y front
set key noautotitle

stats FILE u 0 nooutput   # get number of blocks, i.e. polygons
N = STATS_blocks

getArea(colX,colY)   = ($0==0?(Area=0, x1=column(colX), y1=column(colY)) : 0, \
                        x0=x1, y0=y1, x1=column(colX), y1=column(colY), Area=Area+0.5*(y1+y0)*(x1-x0))
getMinMax(colX,colY) = (x2=column(colX), y2=column(colY), $0==0? (xMin=xMax=x2, yMin=yMax=y2) : \
                       (x2<xMin?xMin=x2:0, x2>xMax?xMax=x2:0, y2<yMin?yMin=y2:0, y2>yMax?yMax=y2:0))
Areas = Centers      = ''
do for [i=1:N] {
    stats FILE u (getArea(1,2),getMinMax(1,2)) index i-1 nooutput
    Areas   = Areas.sprintf(" %g",abs(Area))
    Centers = Centers.sprintf(' %g %g',0.5*(xMin+xMax),0.5*(yMin+yMax))
}
CenterX(n) = real(word(Centers,int(column(n))*2+1))
CenterY(n) = real(word(Centers,int(column(n))*2+2))
Area(n)    = real(word(Areas,int(column(n)+1)))
myColors   = "0xff0000 0x00ff00 0x0000ff"
myColor(i) = sprintf("#%06x",int(word(myColors,(i-1)%words(myColors)+1)))

plot for [i=1:N] FILE u 1:2 index i-1 w filledcurves lc rgb myColor(i), \
     '+' u (CenterX(0)):(CenterY(0)):(sprintf("A=%g",Area(0))) every ::0::N-1 w labels center
### end of script

Result: