Monday, 14 October 2019

Fun with Near-Infrared Drone Imaging

Experimenting with Drone imaging with NIR converted cameras and specialist filters has led me to develop image processing techniques to create NDVI images but it has also allowed me to create visually interesting, (borderline psychedelic!), videos and images.



Travelling with my drone to the beaches of Porto, Portugal allowed me to experience great photo opportunities of the very visually dynamic Atlantic coastline of Portugal.



Choosing different filters also makes sure that you never see the same scene twice, no matter where you are. Even a relatively overcast sky can be very vibrant and beautiful in the Near-Infrared.



Of course the vibrancy of vegetation is the easiest thing to notice in contrast with the relative darkness of other features in an image, such as Indy the cat here!



It is also interesting to notice that the fake grass on the lawn does not reflect strongly in the NIR, unlike the weeds! NIR reflectance is not merely to do with the pigment of vegetation, it is in all likelihood due to the complex structure of the spongy mesophyll layers and the water content in the xylem vessels that run through healthy vegetation, particularly leaves. This is why, as discussed previously, NIR and NDVI is a very valuable indicator of vegetation classification and of plant health.

In any event, whether the application be scientific or artistic, there is always a lot of fun to be had in opening up a new visual perspective on the everyday features of the world.




Notes:


For standard RGB batch processing of DNG (RAW) images files obtained from my DJI Mavic Pro 2 I use the following Python script which uses the rawpy wrapper for the Libraw library for convenience. Here it is:


# -*- coding: utf-8 -*-
"""
Created on Sat Oct 12 03:38:28 2019

@author: cosmi

Working Batch Processer for DNG to JPG

rawpy is an easy-to-use Python wrapper for the LibRaw library.
rawpy works natively with numpy arrays and supports a lot of options,
including direct access to the unprocessed Bayer data
It also contains some extra functionality for finding and repairing hot/dead pixels.
import rawpy.enhance for this


"""

import os
import rawpy
import imageio


for infile in os.listdir("./"):
    print( "file : " + infile)
    if infile[-3:] == "tif" or infile[-3:] == "DNG" :
       # print "is tif or DNG (RAW)"
       outfile = infile[:-3] + "jpg"
       raw = rawpy.imread(infile)
       print( "new filename : " + outfile)
       rgb = raw.postprocess()
       imageio.imsave(outfile, rgb)

Wednesday, 2 October 2019

Specialist NDVI Filter Developments for Mavic Pro 2

I have been developing a brand new series of NDVI specialised filters for use with both standard and IR converted drone cameras. Using these prototypes I am hoping to make some new developments in the use of single-camera near-infrared imaging and NDVI analysis for use in affordable environmental sensing and monitoring worldwide.

The 2 filters are essentially dual bandpass filters, made of high quality glass, with different transmission characteristics in the near-UV to blue band and Near-infrared band.

They are circular filters embedded in a custom housing for use with the Mavic 2 Pro Hasselblad Camera.

Filter#1



Filter#2



With its large sensor area and excellent image vibrancy this sensor can detect the near-infrared band easily for distinguishing plants from other objects in a landscape image.

The following graphs show transmission of the 2 filters.

Filter #1 has the largest transmission in the near-infrared band and has a nice plateau which captures much of the characteristic NIR reflectance pattern from healthy plants.



With IR converted cameras (hot mirror filter removal) the large resulting transitivity with this kind of filter can sometimes overexpose the sensor if taken during high brightness, i.e. during a very sunny day, so in practice although it provides a nice profile the high transmission does not necessarily provide the greatest way to characterize one section of healthy vegetation from another, which is the key factor we want to move forward beyond NIR imaging as just a nice trick for hobbyist photographers and using this technique as a reliable classification tool.

With standard unconverted cameras, (with no hot mirror filter removal) the depleted transmission of the NIR with this kind of filter causes the opposite problem to occur, namely very little IR is recorded by the sensor which means classification is limited. However, with camera settings adjustments, ISO and Exposure settings set to an optimum level for example, we can overcome this and still create NIR profiling so NIR filters can indeed be used by non converted cameras within some limits, i.e. the images must be taken during a very bright clear day at noon to minimize shadowing.


Filter #2 is again a dual bandpass filter but with much narrower and defined transmission lines.



We see the NIR band is far more diminished in transitivity, which does eliminate some of the over-saturation problem with the IR converted cameras. The UV-Blue band is more focused around the near-UV than in filter #1, which means that this filter can pick up more UV information and use it to balance more with the NIR information which could be useful for more novel image analysis beyond the scope of conventional NDVI. These are topics for another article*(see notes)

Image testing is ongoing with these filter designs, with both IR converted and unconverted drone cameras.

In practice I found the characteristics of Filter#2 to be sharper than Images taken with Filter#1. This has advantages in more overcast drone imaging conditions, however in bright daylight conditions Filter#1 tends to be better as shadows become less pronounced by the same criterion.


The following video shows a NIR test flight of a reservoir in Ireland called Blessington Lake, officially Pollaphuca. It is a bird conservation site and source of water for local use. This reservoir was created artificially by the construction of Pollaphuca Dam for use in Hydroelectricity production.



I find this is a very good site to test general environmental boundaries, such as that between a forested area, a shoreline and the water itself.

This test also shows that shadows are diminished if the drone imaging is performed during overcast weather, however the overall NIR reflectance information is very clear. This is an important detail as NDVI generating satellites cannot penetrate cloud covered areas, so Drones add an advantage in this regard.

Moreover, the resolution on the ground using a drone at a height of 100m can be as low as 2cm per pixel, a resolution not possible using satellites and this information is needed for close environmental monitoring and species classification, particularly in wild areas.

Staring off with an image we can begin an example NDVI plot


JPG NIR Image Taken with Mavic 2 Pro, using Filter#1

Either the JPG or DNG (i.e RAW) image format can be converted to NDVI

DNG(RAW) NIR Image Taken with Mavic 2 Pro, using Filter#1


Using the Following Infrablue NDVI Formula:

 

The red and blue colour channels can be used to compute an NDVI colourmap using Matlab (see Appendix for Codes):
JPG NDVI Image (processed in Matlab)

DNG Files contain more information as the colour channels are recorded separately on a 1x3 Matrix and so no splitting is necessary.The resulting image is processed as a PNG or een a TIF image if we like. In any vent we see shadows where NDVI should be zero much more defined as they would have been recorded.

DNG NDVI Image (processed in Matlab)

In the DNG Image the trees are far more defined relative to the grass, the shoreline is also far better defined relative to the water.

Similarly Processed using a standalone Python Script (with a slight amendment in colourmap labeling as to show that the water should indeed be absorbing IR, thus having a negative NDVI):




Updates will be posted here along with a showcasing of videos, images, NDVI batch processing code for use in 3D NDVI Plotting and spectral information on the filters used for ease of method replication.

It is hoped with new filter developments that Enhanced NDVI conversion can be achieved with a single use DJI Mavic Pro 2 Drone Camera, using a filter that also uses the Green Colour Channel for leverage adding to distinguishing different plants, soils, water, algae, rocks and other objects.

I am also hoping that analysis of captured video and making careful notes of exposure settings used in the field in a recent expedition to Portugal will allow me to engineer more refined versions of these filters in the coming weeks.


Appendix:

  • Notes:


Some UV standalone imaging applications, which could be leveraged against simultaneously captured NIR information include imaging fluorescent objects (certain coral, algae, flowering plants, animal species) , recording sulfur aerosols in the air, as sulfur aerosols transmit visible and NIR but block UV. Another possible imaging application is seeing through certain flames (UV and blue light travels through orange portion of flames and this has been used by the steel industry in welding for decades as a proven technology to image an object that is shrouded by embers or sparks). Some of these imaging techniques, within reason, could be indexed in a similar way to NDVI. 


  • Matlab Code#1 (JPG Input Only):

%Experimental Vegetation Index Mapping program using DJI Mavic Pro DNG 16-bit images taken using InfraBlue Filter 
%Works with JPG Images only - DNG Files need a separate script
%(c)-J. Campbell MuonRay Enterprises Drone-based Vegetation Index Project 2017-2019 
% Read original image.
rgbImage = imread('DJI_0182.JPG');
% Extract colour channels.
redChannel = rgbImage(:,:,1); % Red channel
greenChannel = rgbImage(:,:,2); % Green channel
blueChannel = rgbImage(:,:,3); % Blue channel

figure(1);

subplot(2,2,1)
imshow(redChannel)
title('Red Channel')

subplot(2,2,2)
imshow(greenChannel)
title('Green Channel')


subplot(2,2,3)
imshow(blueChannel)
title('Blue Channel')

B = im2single(blueChannel); %VIS Channel Blue - which is used as a magnitude scale 
R = im2single(redChannel); %NIR Chaneel RED which is the IR reflectance of healthy vegeation 


InfraBlue = (R - B)./(R + B);
InfraBlue = double(InfraBlue);



%% Stretch NDVI to 0-255 and convert to 8-bit unsigned integer
InfraBlue = floor((InfraBlue + 1) * 32767); % [-1 1] -> [0 256] for 8-bit display range(*128), [0, 65535] for 16-bit display range (*32500)
InfraBlue(InfraBlue < 0) = 0;             % not really necessary, just in case & for symmetry
InfraBlue(InfraBlue > 65535) = 65535;         % in case the original value was exactly 1
InfraBlue = uint16(round(InfraBlue));             % change data type from double to uint16
% InfraBlue = uint8(InfraBlue);             % change data type from double to uint8

NDVImag = double(InfraBlue);


figure(3);

myColorMap = jet(65535); % Whatever you want.
rgbImage = ind2rgb(NDVImag, myColorMap);

imagesc(rgbImage,[0 1]), colorbar

%imshow(rgbImage), colorbar
%c = colorbar % Add a color bar to indicate the scaling of color


imwrite(rgbImage, 'InfraBlueJPGInput.png');


  • Matlab Code#2 (DNG Input Only):
  • Main File:


%Experimental Vegetation Index Mapping program using DJI Mavic Pro DNG 16-bit images taken using InfraBlue Filter 
%Works with DNG Images only - needs loadDNG script in same directory
%(c)-J. Campbell MuonRay Enterprises Drone-based Vegetation Index Project 2017-2019 

rawData = loadDNG('DJI_0182.DNG'); % load it "functionally" from the command line

%DNG is just a fancy form of TIFF. So we can read DNG with TIFF class in MATLAB

warning off MATLAB:tifflib:TIFFReadDirectory:libraryWarning


%Transpose Matrices for image data because DNG image format is row-major, and Matlab matrix format is column-major.

%rawData = permute(rawData, [2 1 3]);

%Assume Bayer mosaic sensor alignment.
%Seperate to mosaic components.
J1 = rawData(1:2:end, 1:2:end);
J2 = rawData(1:2:end, 2:2:end);
J3 = rawData(2:2:end, 1:2:end);
J4 = rawData(2:2:end, 2:2:end);


figure(1);

subplot(2,2,1)
imshow(J1)
title('Red Channel')

subplot(2,2,2)
imshow(J2)
title('Green Channel 1')

subplot(2,2,3)
imshow(J3)
title('Green Channel 2')

subplot(2,2,4)
imshow(J4)
title('Blue Channel')


J1 = imadjust(J1, [0.09, 0.91], [], 0.45); %Adjust image intensity, and use gamma value of 0.45

J1(:,:,1) = J1(:,:,1)*0.80; %Reduce the overall colour temperature by factor of 0.80

B = im2single(J4(:,:,1)); %VIS Channel Blue - which is used as a magnitude scale 
R = im2single(J1(:,:,1)); %NIR Chaneel RED which is the IR reflectance of healthy vegeation 


InfraBlue = (R - B)./(R + B);
InfraBlue = double(InfraBlue);


%% Stretch NDVI to 0-255 and convert to 8-bit unsigned integer
InfraBlue = floor((InfraBlue + 1) * 32767); % [-1 1] -> [0 256] for 8-bit display range(*128), [0, 65535] for 16-bit display range (*32500)
InfraBlue(InfraBlue < 0) = 0;             % not really necessary, just in case & for symmetry
InfraBlue(InfraBlue > 65535) = 65535;         % in case the original value was exactly 1
InfraBlue = uint16(round(InfraBlue));             % change data type from double to uint16
% InfraBlue = uint8(InfraBlue);             % change data type from double to uint8

NDVImag = double(InfraBlue);


figure(3);

myColorMap = jet(65535); % Whatever you want.
rgbImage = ind2rgb(NDVImag, myColorMap);

imagesc(rgbImage,[0 1]), colorbar

%imshow(rgbImage), colorbar
%c = colorbar % Add a color bar to indicate the scaling of color


imwrite(rgbImage, 'InfraBlueFin.png');

  • loadDNG Header File:

function [rawData, tinfo]= loadDNG(dngFilename)   
    if(exist(dngFilename,'file'))
        tinfo = imfinfo(dngFilename);
        t = Tiff(dngFilename,'r');
        offsets = getTag(t,'SubIFD');
        setSubDirectory(t,offsets(1));
        rawData = t.read();
        t.close();
    else
        if(nargin<1 || isempty(dngFilename))
            dngFilename = 'File';
        end
        fprintf(1,'%s could not be found\n',dngFilename);
        rawData = [];
    end
end



  • Python Code (JPG Input Only)


"""
Experimental Vegetation Index Mapping program using DJI Mavic 2 Pro
JPEG 16-bit combo images taken using InfraBlue Filter

%(c)-J. Campbell MuonRay Enterprises 2019
This Python script was created using the Spyder Editor

"""

import warnings
warnings.filterwarnings('ignore')
from scipy import misc
import imageio
import numpy as np
from matplotlib import pyplot as plt  # For image viewing

#!/usr/bin/python
import getopt
import sys

import matplotlib.pyplot as plt
from matplotlib import colors
from matplotlib import ticker
from matplotlib.colors import LinearSegmentedColormap

#dng reading requires libraw to work

# Open an image
image = misc.imread('DJI_0182.JPG')

# Get the red band from the rgb image, and open it as a numpy matrix
#NIR = image[:, :, 0]
         
#ir = np.asarray(NIR, float)


ir = (image[:,:,0]).astype('float')


# Get one of the IR image bands (all bands should be same)
#blue = image[:, :, 2]

#r = np.asarray(blue, float)

r = (image[:,:,2]).astype('float')


# Create a numpy matrix of zeros to hold the calculated NDVI values for each pixel
ndvi = np.zeros(r.size)  # The NDVI image will be the same size as the input image

# Calculate NDVI
ndvi = np.true_divide(np.subtract(ir, r), np.add(ir, r))


# Display the results
output_name = 'InfraBlueNDVI3.jpg'

#a nice selection of grayscale colour palettes
cols1 = ['blue', 'green', 'yellow', 'red']
cols2 =  ['gray', 'gray', 'red', 'yellow', 'green']
cols3 = ['gray', 'blue', 'green', 'yellow', 'red']

cols4 = ['black', 'gray', 'blue', 'green', 'yellow', 'red']

def create_colormap(args):
    return LinearSegmentedColormap.from_list(name='custom1', colors=cols3)

#colour bar to match grayscale units
def create_colorbar(fig, image):
        position = fig.add_axes([0.125, 0.19, 0.2, 0.05])
        norm = colors.Normalize(vmin=-1., vmax=1.)
        cbar = plt.colorbar(image,
                            cax=position,
                            orientation='horizontal',
                            norm=norm)
        cbar.ax.tick_params(labelsize=6)
        tick_locator = ticker.MaxNLocator(nbins=3)
        cbar.locator = tick_locator
        cbar.update_ticks()
        cbar.set_label("NDVI", fontsize=10, x=0.5, y=0.5, labelpad=-25)

fig, ax = plt.subplots()
image = ax.imshow(ndvi, cmap=create_colormap(colors))
plt.axis('off')

create_colorbar(fig, image)

extent = ax.get_window_extent().transformed(fig.dpi_scale_trans.inverted())
fig.savefig(output_name, dpi=600, transparent=True, bbox_inches=extent, pad_inches=0)
        # plt.show()