In this post, I will face a difficult challenge: detect whether a plant is flourishing or not. At this moment in time, this feature requires some more experiments and works to make a reliable solution, but I think the solution is promising
NDVI
Vegetation is green because plant leaves reflect green light -- they don't use much of it for photosynthesis (see Figure 1).
Instead they use lots of the blue and red wavelengths in sunlight. The pigments in leaves absorb this light to power photosynthesis which converts CO2, water, and nutrients into carbohydrates (food). In general, you can estimate the productivity or vigor of vegetation by how much blue and red light it is absorbing. Photosynthetic pigments do not use the longer, invisible wavelengths of near infrared light and reflect almost all of it away (this helps prevent the leaves from overheating).
Shortly after the launch of the first Landsat satellite in 1972, scientists began using the data from its sensors to estimate the productivity of vegetation by comparing the amount of red light reflected (there is not much from healthy plants) to the amount of near infrared light reflected (there is a lot). The amount of infrared light reflected from vegetation is a good indicator of how bright the sunlight was at any moment (leaves reflect almost all IR). Comparing that to the amount of reflected red light can suggest what proportion of the sunlight was being absorbed by the plants.
That relationship is a good measure of the amount of photosynthetically active biomass. They quickly settled on an index of plant productivity called NDVI for Normalized Difference Vegetation Index defined as
N - R / N + R
where N is the amount of radiation in the near-infrared spectrum and R is the amount of radiation in the red spectrum.
It is possible to capture all the information needed to compute NDVI in just one camera using a so-called "superblue" filter. Spectral transmission of a "superblue" filter is shown in the image below. Most red light (600 to 700 nm) is blocked, but other visible light and most near infrared light (> 700 nm) passes through. This type of filter can replace the IR block filter to allow NDVI to be produced from a single camera. The red channel in the image file will record mostly infrared light, and the blue channel will record light that plants use for photosynthesis.
After some theory, nothing is better than some coding...
NDVI Implementation
First of all, I need to split the original RGB image (actually an NGB, because R channel is replaces by the Near-infrared radiation) into its three components
CvSize imgSize = cvGetSize(image);
IplImage* n = cvCreateImage(imgSize, image->depth, 1);
IplImage* g = cvCreateImage(imgSize, image->depth, 1);
IplImage* b = cvCreateImage(imgSize, image->depth, 1);
cvSplit(image, b, g, n, NULL);
Then I created the image where the NDVI will be shown. This has only 1 channel
IplImage* imgNDVI = cvCreateImage(imgSize, image->depth, 1);
I will also map the NDVI range of values (from -1 to 1) to RGB colors. The resulting image will have three channels.
IplImage* imgOut = cvCreateImage(imgSize, image->depth, 3);
Next I need some pointer arithmetic to make image pixels access faster
unsigned int stepIn = n->widthStep / (n->depth / 8);
unsigned int stepOut = imgOut->widthStep / (imgOut->depth / 8);
unsigned int imgIn_width = n->width;
unsigned int imgIn_height = n->height;
unsigned int imgIn_offset = 0;
unsigned int imgOut_width = imgOut->width;
unsigned int imgOut_height = imgOut->height;
unsigned int imgOut_offset = 0;
if(n->roi)
{
imgIn_width = n->roi->width;
imgIn_height = n->roi->height;
imgIn_offset = n->roi->xOffset + (n->roi->yOffset * stepIn);
}
if(imgOut->roi)
{
imgOut_width = imgOut->roi->width;
imgOut_height = imgOut->roi->height;
imgOut_offset = imgOut->roi->xOffset + (imgOut->roi->yOffset * stepOut);
}
unsigned char *imgDataN = (unsigned char *)n->imageData + imgIn_offset;
unsigned char *imgDataB = (unsigned char *)b->imageData + imgIn_offset;
unsigned char *imgDataOut = (unsigned char *)imgOut->imageData + imgOut_offset;
unsigned char *imgDataNDVI = (unsigned char *)imgNDVI->imageData + imgIn_offset;
#define imageN(X, Y) imgDataN[(X) + (Y)*stepIn]
#define imageB(X, Y) imgDataB[(X) + (Y)*stepIn]
#define imageOut(X, Y, ch) imgDataOut[(X)*3 + (Y)*stepOut + ch]
#define imageNDVI(X, Y) imgDataNDVI[(X) + (Y)*stepIn]
All the above stuff is just to compute the number of bytes for row in each image involved in the NDVI calculation. The three macro imageN, imageB and imageOut make it easy to access the image pixels. Thanks to this macros, I can write code like
imageN(10, 100)
to access the pixel value at position (10, 100). Note that the imageOut has an extra parameter (ch) which is the channel. Because OpenCV uses the BGR pixel format, ch's permitted values are 0 (for B channel), 1 (for G channel) and 2 (for R channel)
To compute NDVI, I just loop on each pixel and convert the NDVI value to the corresponding RGB color
#define NVDI_MIN -1
#define NVDI_MAX +1
#define PIXEL_MIN 0
#define PIXEL_MAX 255
// for each pixel, compute (N-B) / (N+B)
for (unsigned int y=0; y<imgIn_height; y++)
{
for (unsigned int x=0; x<imgIn_width; x++)
{
double N = imageN(x, y);
double B = imageB(x, y);
double NVDI = (N-B) / (N+B);
double px = (NVDI - NVDI_MIN) * (PIXEL_MAX-PIXEL_MIN) / (NVDI_MAX-NVDI_MIN) + PIXEL_MIN;
imageNDVI(x,y) = (unsigned char)px;
uint8_t r, g, b;
ValueToColor(NVDI, -1, +1, &r, &g, &b);
imageOut(x,y, 0) = b;
imageOut(x,y, 1) = g;
imageOut(x,y, 2) = r;
}
}
The imgNDVI is then searched for blobs that represents flourishing plants
IplImage* frame = cvCreateImage(imgSize, imgOut->depth, imgOut->nChannels);
cvConvertScale(imgOut, frame, 1, 0);
double threshold = (0.05 - NVDI_MIN) * (PIXEL_MAX-PIXEL_MIN) / (NVDI_MAX-NVDI_MIN) + PIXEL_MIN;
IplImage* segmentated = cvCreateImage(imgSize, 8, 1);
cvInRangeS(imgNDVI, CV_RGB((uint8_t)threshold, (uint8_t)threshold, (uint8_t)threshold),
CV_RGB(255, 255, 255), segmentated);
IplImage* labelImg = cvCreateImage(cvGetSize(image), IPL_DEPTH_LABEL, 1);
CvBlobs blobs;
cvLabel(segmentated, labelImg, blobs);
cvFilterByArea(blobs, 200, 1000000);
result = blobs.size();
Results and videos in my next post...
