Changed Composite back to Original

Canadian-Severe-Storms-Laboratory · Oct 2, 2024 · 8d10e9f · 8d10e9f
1 parent 5266c0f
commit 8d10e9f
Show file tree

Hide file tree

Showing 6 changed files with 147 additions and 96 deletions.
diff --git a/LandCover.cs b/LandCover.cs
@@ -7,16 +7,19 @@
 using ArcGIS.Desktop.Framework.Threading.Tasks;
 using System.Threading.Tasks;
 using System.Collections.Generic;
+using System.IO.Packaging;
 
 public static class LandCover
 {
 
     #if DEBUG
-        private static string path = @"C:\Users\danie\Documents\Experiments\SAR analysis\landcover-2020-classification.tif";
+        private static string path = @"C:\Users\danie\Documents\Experiments\SAR analysis";
     #else
-        private static string path = AddinAssemblyLocation() + "\\landcover-2020-classification.tif";
+        private static string path = AddinAssemblyLocation();
     #endif
 
+    private static string fileName = "landcover-2020-classification.tif";
+
     public enum LandCoverType
     {
         ConiferousForest = 1,
@@ -40,7 +43,7 @@ public enum LandCoverType
 
     public static Mat GetSection(RasterLayer referenceLayer) {
 
-        RasterLayer layer = LoadRasterLayer(path);
+        RasterLayer layer = LoadRasterLayer(path, fileName);
 
         Raster raster = layer.GetRaster();
         Raster refRaster = referenceLayer.GetRaster();

diff --git a/Main.xaml b/Main.xaml
@@ -41,8 +41,11 @@
             </TabItem>
             <TabItem Header="Results" Width="100">
                 <Grid>
-                    <WpfPlot x:Name="resultsPlot" Margin="10,10,10,10"/>
-
+                    <CheckBox x:Name="histogramToggle" Style="{DynamicResource Esri_CheckboxToggleSwitch}" HorizontalAlignment="Right" VerticalAlignment="Top" Margin="0,8,50,0" Unchecked="histogramToggle_Unchecked" Checked="histogramToggle_Checked"/>
+                    <Label Content="VH" HorizontalAlignment="Right" Margin="0,1,80,0" VerticalAlignment="Top"/>
+                    <Label Content="VV" HorizontalAlignment="Right" Margin="0,1,30,0" VerticalAlignment="Top"/>
+                    <WpfPlot x:Name="VHPlot" Margin="10,30,10,10"/>
+                    <WpfPlot x:Name="VVPlot" Margin="10,30,10,10" Visibility="Hidden"/>
                 </Grid>
             </TabItem>
             <TabItem Header="Console" Width="100">

diff --git a/Main.xaml.cs b/Main.xaml.cs
@@ -1,6 +1,7 @@
 using ArcGIS.Core.Geometry;
 using ArcGIS.Desktop.Framework.Threading.Tasks;
 using ArcGIS.Desktop.Mapping;
+using ScottPlot;
 using System;
 using System.Collections.Generic;
 using System.Threading.Tasks;
@@ -36,32 +37,20 @@ private async void RunAnalysis(object sender, RoutedEventArgs e)
 
                 SarAnalyser sarAnalyser = BuildSarAnalyser();
 
-                List<double> values = new();
+                List<double> VHvalues = new();
+                List<double> VVvalues = new();
 
                 await QueuedTask.Run(async () =>
                 {
                     PrintTitle();
 
-                    values = await sarAnalyser.Analyze();
+                    (VHvalues, VVvalues) = await sarAnalyser.Analyze();
 
                     Console.WriteLine("Done.\n");
                 });
 
-                var plot = resultsPlot.Plot;
-                ScottPlot.Statistics.BasicStats stats = new ScottPlot.Statistics.BasicStats(values.ToArray());
-                ScottPlot.Statistics.Histogram hist = new(min: stats.Min, stats.Max, binCount: 10);
-                hist.AddRange(values);
-
-                var barPlot = plot.AddBar(values: hist.GetProbability(), positions: hist.Bins);
-                barPlot.BarWidth = (stats.Max - stats.Min) / 10;
-                plot.YAxis.Label("Density");
-                plot.XAxis.Label("Pixel Value");
-                plot.Title("Normalized VH*VV difference histogram for AOI");
-                plot.SetAxisLimits(yMin: 0);
-                plot.AddVerticalLine(x: stats.Mean, color: System.Drawing.Color.Black, style: ScottPlot.LineStyle.Dash, label: "Mean: " + Math.Round(stats.Mean, 2).ToString());
-                plot.Legend(location: ScottPlot.Alignment.UpperRight);
-
-                resultsPlot.Refresh();
+                PlotHistogram(VHvalues, VHPlot, "VH");
+                PlotHistogram(VVvalues, VVPlot, "VV");
 
                 tabController.SelectedIndex = 2;
 
@@ -76,6 +65,29 @@ await QueuedTask.Run(async () =>
             runButton.IsEnabled = true;
         }
 
+        private void PlotHistogram(List<double> values, WpfPlot plot, string bandName)
+        {
+            values.Sort();
+            double min = values[(int)(values.Count * 0.01)];
+            double max = values[(int)(values.Count * 0.99)];
+
+            var plt = plot.Plot;
+            ScottPlot.Statistics.BasicStats stats = new ScottPlot.Statistics.BasicStats(values.ToArray());
+            ScottPlot.Statistics.Histogram hist = new(min, max, binCount: 8, addOutliersToEdgeBins: true);
+            hist.AddRange(values);
+
+            var barPlot = plt.AddBar(values: hist.GetProbability(), positions: hist.Bins);
+            barPlot.BarWidth = (max - min) / 8;
+            plt.YAxis.Label("Percentage of Pixels");
+            plt.XAxis.Label("Pixel Value");
+            plt.Title(bandName + " difference histogram for AOI");
+            plt.SetAxisLimits(yMin: 0);
+            plt.AddVerticalLine(x: stats.Mean, color: System.Drawing.Color.Red, width: 2, style: ScottPlot.LineStyle.Dash, label: "Mean: " + Math.Round(stats.Mean, 2).ToString());
+            plt.Legend(location: ScottPlot.Alignment.UpperRight);
+
+            plot.Refresh();
+        }
+
         private SarAnalyser BuildSarAnalyser()
         {
             return new SarAnalyser(preEventSelection.GetFilePath(), postEventSelection.GetFilePath(), aoiMaskSelection.GetSelectedLayer());
@@ -124,5 +136,17 @@ private void PrintTitle()
                           "/_/  \\___/_/ /_//_/\\_,_/\\_,_/\\___/ /___/_/ |_/_/|_| \n\n" +
                           "----------------------------------------------------------\n\n");
         }
+
+        private void histogramToggle_Unchecked(object sender, RoutedEventArgs e)
+        {
+            VHPlot.Visibility = Visibility.Visible;
+            VVPlot.Visibility = Visibility.Hidden;
+        }
+
+        private void histogramToggle_Checked(object sender, RoutedEventArgs e)
+        {
+            VVPlot.Visibility = Visibility.Visible;
+            VHPlot.Visibility = Visibility.Hidden;
+        }
     }
 }
diff --git a/SarAnalyser.cs b/SarAnalyser.cs
@@ -1,18 +1,12 @@
 using System;
 using System.Collections.Generic;
-using System.Linq;
-using System.Text;
 using System.Threading.Tasks;
-using System.Windows;
 using ArcGIS.Core.CIM;
-using ArcGIS.Core.Data;
+using ArcGIS.Core.Data.Raster;
 using ArcGIS.Core.Geometry;
-using ArcGIS.Desktop.Core.Geoprocessing;
-using ArcGIS.Desktop.Framework.Threading.Tasks;
+using ArcGIS.Desktop.Internal.DesktopService;
 using ArcGIS.Desktop.Mapping;
 using OpenCvSharp;
-using ScottPlot;
-using static System.Runtime.InteropServices.JavaScript.JSType;
 using static ArcGISUtils.Utils;
 using static LandCover.LandCoverType;
 
@@ -39,16 +33,20 @@ public SarAnalyser(string prePath, string postPath, FeatureLayer aoiLayer)
             System.IO.Directory.CreateDirectory(resultsPath);
         }
 
-        public async Task<List<double>> Analyze()
+        public async Task<Tuple<List<double>, List<double>>> Analyze()
         {
-            var (collocatedLayer, compositeLayer) = ContructSARRasters();
+            var collocatedLayer = ContructSARRasters();
 
+            Console.WriteLine("Loading Collocated Raster");
             List<Mat> bands = LoadCollocatedRaster(collocatedLayer);
 
-            Mat diffMul = Normalized_VH_VV_Difference(bands);
+            Console.WriteLine("Constructing Composite Raster");
+            Mat diff = ReverseDifference(bands);
 
-            WriteComposite(diffMul, compositeLayer);
+            Console.WriteLine("Loading Composite Raster");
+            RasterLayer compositeLayer = LoadComposite(diff);
 
+            Console.WriteLine("Computing AOI Statistics");
             Mat mask = PolygonAOIMask(collocatedLayer);
 
             //Mat LandCoverMask = LandCover.TypeMask(LandCover.GetSection(collocatedLayer), [ConiferousForest, TaigaForest, DeciduousForest, MixedForest, Grassland, Shrubland]);
@@ -57,20 +55,42 @@ public async Task<List<double>> Analyze()
 
             //WriteComposite(LandCoverMask.Resize(polyMask.Size(), interpolation: InterpolationFlags.Nearest), compositeLayer);
 
+            return GetAOIValues(mask, diff);
+        }
+
+
+        private Tuple<List<double>, List<double>> GetAOIValues(Mat mask, Mat diff)
+        {
             mask.GetArray(out byte[] maskArray);
-            diffMul.GetArray(out float[] pixels);
+            Cv2.Split(diff, out Mat[] diffs);
+            diffs[0].GetArray(out float[] VHpixels);
+            diffs[1].GetArray(out float[] VVpixels);
 
-            List<double> values = new();
+            List<double> VHvalues = new(VHpixels.Length);
+            List<double> VVvalues = new(VVpixels.Length);
 
             for (int i = 0; i < maskArray.Length; i++)
             {
-                if (maskArray[i] > 0 && pixels[i] > 0)
+                if (maskArray[i] > 0)
                 {
-                    values.Add(pixels[i]);
+                    VHvalues.Add(VHpixels[i]);
+                    VVvalues.Add(VVpixels[i]);
                 }
             }
 
-            return values;
+            return new (VHvalues, VVvalues);
+        }
+
+        private Mat ReverseDifference(List<Mat> bands)
+        {
+            Mat diffVH = bands[0].Subtract(bands[2]);
+            Mat diffVV = bands[1].Subtract(bands[3]);
+
+            Mat diff = new Mat();
+
+            Cv2.Merge([diffVH, diffVV], diff);
+
+            return diff;
         }
 
         private Mat PolygonAOIMask(RasterLayer collocatedLayer)
@@ -96,18 +116,25 @@ private string getStandardAOIExtentString()
             return "\"POLYGON((" + xmin + " " + ymax + ", " + xmax + " " + ymax + ", " + xmax + " " + ymin + ", " + xmin + " " + ymin + ", " + xmin + " " + ymax + "))\"";
         }
 
-        private void WriteComposite(Mat diff, RasterLayer compositeLayer)
+        private RasterLayer LoadComposite(Mat diff)
         {
-            WriteRaster<float>(diff, compositeLayer.GetRaster());
+            var compositeDataset = OpenRasterDataset(resultsPath, "composite.tif");
+            Raster compositeRaster = compositeDataset.CreateFullRaster();
 
-            var colorizer = compositeLayer.GetColorizer() as CIMRasterRGBColorizer;
-            colorizer.RedBandIndex = 1;
-            colorizer.GreenBandIndex = 2;
+            WriteRaster<float>(diff, compositeRaster);
+
+            var compositeLayer = LoadRasterLayer(resultsPath, "composite.tif");
+
+            CIMRasterRGBColorizer colorizer = new();
+            colorizer.RedBandIndex = 2;
+            colorizer.GreenBandIndex = 3;
             colorizer.BlueBandIndex = 0;
-            colorizer.StretchType = RasterStretchType.MinimumMaximum;
+            colorizer.DisplayBackgroundValue = true;
             compositeLayer.SetColorizer(colorizer);
 
             MapView.Active.Redraw(true);
+
+            return compositeLayer;
         }
 
         private List<Mat> LoadCollocatedRaster(RasterLayer collocatedLayer)
@@ -128,6 +155,7 @@ private List<Mat> LoadCollocatedRaster(RasterLayer collocatedLayer)
             colorizer.GreenBandIndex = 2;
             colorizer.BlueBandIndex = 3;
             colorizer.DisplayBackgroundValue = true;
+
             collocatedLayer.SetColorizer(colorizer);
 
             List<Mat> bands = [postBands[0], postBands[1], preBands[0], preBands[1]];
@@ -140,22 +168,23 @@ private List<Mat> LoadCollocatedRaster(RasterLayer collocatedLayer)
             return bands;
         }
 
-        private Tuple<RasterLayer, RasterLayer> ContructSARRasters()
+        private RasterLayer ContructSARRasters()
         {
+            Console.WriteLine("Constructing Pre Raster...");
             string preCorrectedPath = CorrectSARData(prePath, "pre");
+            Console.WriteLine("\nConstructing Post Raster...");
             string postCorrectedPath = CorrectSARData(postPath, "post");
 
             string collocatedPath = "\"" + resultsPath + "\\collocated.dim\"";
 
+            Console.WriteLine("\nConstructing Collocated Raster...");
             SnapWrapper.RunCommand("collocate " + postCorrectedPath + " " + preCorrectedPath + " -t " + collocatedPath);
             SnapWrapper.RunCommand("speckle-filter " + collocatedPath + " -Pfilter=\"Refined Lee\" -PsourceBands=\"Sigma0_VH_M, Sigma0_VV_M, Sigma0_VH_S, Sigma0_VV_S\" -f Geotiff -t \"" + resultsPath + "\\collocated.tif\"");
-            SnapWrapper.RunCommand("speckle-filter " + collocatedPath + " -Pfilter=\"Refined Lee\" -PsourceBands=\"Sigma0_VH_M, Sigma0_VV_M\" -f Geotiff -t \"" + resultsPath + "\\composite.tif\"");
-
-            RasterLayer collocatedLayer = LoadRasterLayer(resultsPath + "\\collocated.tif");
+            SnapWrapper.RunCommand("speckle-filter " + collocatedPath + " -Pfilter=\"Refined Lee\" -PsourceBands=\"Sigma0_VH_S, Sigma0_VH_M, Sigma0_VV_M\" -f Geotiff -t \"" + resultsPath + "\\composite.tif\"", verbose: false);
 
-            RasterLayer compositeLayer = LoadRasterLayer(resultsPath + "\\composite.tif");
+            RasterLayer collocatedLayer = LoadRasterLayer(resultsPath, "collocated.tif");
 
-            return new(collocatedLayer, compositeLayer);
+            return collocatedLayer;
         }
 
         private string CorrectSARData(string sarPath, string name)

diff --git a/SnapWrapper.cs b/SnapWrapper.cs
@@ -3,9 +3,10 @@
 
 public static class SnapWrapper
 {
-    public static void RunCommand(string command)
+    public static void RunCommand(string command, bool verbose=true)
     {
-        Console.Write("Running SNAP " + command.Split(' ')[0]);
+        if (verbose) Console.Write("   Running SNAP " + command.Split(' ')[0]);
+
         ProcessStartInfo ProcessInfo = new ProcessStartInfo("cmd.exe", "/c gpt " + command);
         ProcessInfo.CreateNoWindow = true;
         ProcessInfo.UseShellExecute = false;