AE中栅格计算的问题

http://www.gisschool.com.tw/Discuz7/redirect.php?tid=119&goto=lastpost

話說今天有人問我關於影像 Raster 計算問題..
給偶了一串很長的公式, 問我要怎麼方便計算..
像 10 - ( ALayer * BLayer ) / CLayer + DLayer * 0.11 + tan ( ELayer ) * cos ( FLayer )
偶就跟他說用 Spatial Analysis 中的 Raster Calculator  的方式去計算就 ok 了ㄚ.

可是對方跟偶說是要動態的, 要去設計視窗給 user 作. 圖層不一定..
好吧.. 就想說之前有作過 透過 MathOp 的方式來設計這串公式.
之前偷是透過 MathOp 這個 function 搞 Times 乘 Plus 加, 當然也有 Tan , Cos 等三角函數...
由於 MathOp 的參數都是兩個 GeoDataSet 的物件, 針對於 10 參數的話..
一種就是建一個 const 為 10 的 Raster  資料.. 再去加減..
或是透過 IPixelBlock 的方式去計算影像的資料.. 
搞的規模很龐大.. 搞的天昏地案日月無光. 搞的一個頭兩個大...
程式就越寫越多. function 就越搞越雜..
寫個漏漏長的幾千行程式後發現.... 發現....

底下是很笨很笨的下去慢慢算..
        Dim pMathOp As IMathOp
        Set pMathOp = New RasterMathOps

        Dim pResultDataset As IGeoDataset
        Set pResultDataset = New RasterDataset
        Set pResultDataset = pMathOp.Times(pRMLayer.Raster, pKMLayer.Raster)
        Set pResultDataset = pMathOp.Div(pResultDataset, pSLayer.Raster)
        Set pResultDataset = pMathOp.ATan(pResultDataset)

然後針對那些參數的加減乘除..
Public Function RasterParamCalc(pResultDataset As IGeoDataset, ctype As String, fnum As Double) As IGeoDataset
        
         Dim pTmpLayer As IRasterLayer
         Set pTmpLayer = New RasterLayer
         pTmpLayer.CreateFromRaster pResultDataset
         
         Dim pRasProps As IRasterProps
         Set pRasProps = pTmpLayer.Raster

         Dim pPB As IPixelBlock3
         Dim pPnt As IPnt
         Set pPnt = New Pnt
         pPnt.SetCoords pRasProps.Width, pRasProps.Height
         Set pPB = pTmpLayer.Raster.CreatePixelBlock(pPnt)
         
         If pPB Is Nothing Then
            MsgBox "影像資料無法取得PixelBlock.."
            Exit Function
         End If
         
         'Dim pRawPixels As IRawPixels
         pPnt.X = 0
         pPnt.Y = 0
         pTmpLayer.Raster.Read pPnt, pPB
                       
             Dim v As Variant
             v = pPB.PixelData(0)
             Dim m, n As Integer
             For m = 0 To pPB.Width - 1
                For n = 0 To pPB.Height - 1
                    If (v(m, n) < 65535 And v(m, n) > -65535) Then
                      If ctype = "+" Then Set v(m, n) = v(m, n) + CDbl(fnum)
                      If ctype = "-" Then Set v(m, n) = v(m, n) - CDbl(fnum)
                      If ctype = "*" Then v(m, n) = v(m, n) * CDbl(fnum)
                      If ctype = "/" Then v(m, n) = v(m, n) / CDbl(fnum)
                    End If
                Next n
             Next m
             pPB.PixelData(0) = v
             Set v = Nothing

        'Write the pixel block with a (0,0) offset of the upper left corner
         Dim pRasterEdit As IRasterEdit
         Set pRasterEdit = pTmpLayer.Raster
         
         If pRasterEdit Is Nothing Then
            MsgBox "影像資料無法取得編輯狀態.."
            Exit Function
         End If
         
         pPnt.X = 0
         pPnt.Y = 0
         pRasterEdit.Write pPnt, pPB
         
         Set RasterParamCalc = pRasterEdit

         Set pTmpLayer = Nothing
         Set pRasProps = Nothing
         Set pPB = Nothing
         Set pPnt = Nothing
         Set v = Nothing
         Set pRasterEdit = Nothing

End Function

忙了好好好好長的一陣子後發現.. 有個神奇的 IMapAlgebraOp
就像 RatserCalculate  可以直接就給他算下企...
天ㄚ.. 我是在作啥...

        Dim pMapAlgebraOp As IMapAlgebraOp
        Set pMapAlgebraOp = New RasterMapAlgebraOp

        pMapAlgebraOp.BindRaster pALayer, "A"
        pMapAlgebraOp.BindRaster pBLayer, "B"
        pMapAlgebraOp.BindRaster pCLayer, "C"
        pMapAlgebraOp.BindRaster pELayer, "E"        
        pMapAlgebraOp.BindRaster pFLayer, "F"

        Dim pCalcDS As IGeoDataset
        Set pCalcDS = New RasterDataset
                
        Set pCalcDS = pMapAlgebraOp.Execute(
                      " 1 + ( [A] / ( ( [E] + 10 ) * [I] * [D] * Cos([F]) * ATan( [B])) ")

C#代码：

private void RasterCaculator(string path) 

        {

            IMap map = MainMap.ActiveView.FocusMap;

            ILayer layer1 = map.get_Layer(0);

            ILayer layer2 = map.get_Layer(1);

            ILayer layer3 = map.get_Layer(2);

            ILayer layer4 = map.get_Layer(3);

            ILayer layer5 = map.get_Layer(4);

            IRasterLayer rasterlayer1 = layer1 as IRasterLayer;

            IRasterLayer rasterlayer2 = layer2 as IRasterLayer;

            IRasterLayer rasterlayer3 = layer3 as IRasterLayer;

            IRasterLayer rasterlayer4 = layer4 as IRasterLayer;

            IRasterLayer rasterlayer5 = layer5 as IRasterLayer;

            IRaster pRaster1 = rasterlayer1.Raster;//获得已知栅格图层pRasterLyr的栅格对象

            IRaster pRaster2 = rasterlayer2.Raster;

            IRaster pRaster3 = rasterlayer3.Raster;

            IRaster pRaster4 = rasterlayer4.Raster;

            IRaster pRaster5 = rasterlayer5.Raster;

            IGeoDataset tempGeodata1 = pRaster1 as IGeoDataset;//

            IGeoDataset tempGeodata2 = pRaster2 as IGeoDataset;

            IGeoDataset tempGeodata3 = pRaster3 as IGeoDataset;

            IGeoDataset tempGeodata4 = pRaster4 as IGeoDataset;

            IGeoDataset tempGeodata5 = pRaster5 as IGeoDataset;

            IMapAlgebraOp rsalgebra = new RasterMapAlgebraOpClass();

            //设置栅格运算空间

            IRasterAnalysisEnvironment rasAnaEnv = (IRasterAnalysisEnvironment)rsalgebra;

            IWorkspaceFactory wsf = new RasterWorkspaceFactoryClass();

            IWorkspace ws=wsf.OpenFromFile("",0);//设置输出空间

            rasAnaEnv.OutWorkspace = ws;

            rsalgebra.BindRaster(tempGeodata1, "1");

            rsalgebra.BindRaster(tempGeodata2, "2");

            rsalgebra.BindRaster(tempGeodata3, "3");

            rsalgebra.BindRaster(tempGeodata4, "4");

            rsalgebra.BindRaster(tempGeodata5, "5");

            IGeoDataset outGetDataset=rsalgebra.Execute("[1] * [2] * [3] * [4] * [5]");

            IRasterLayer pCreatRalyr=new RasterLayerClass();

            pCreatRalyr.CreateFromRaster((IRaster)outGetDataset);

            MainMap.AddLayer(pCreatRalyr);

}