GMT三維繪圖有BUG? 修復它！

阿新 • • 發佈：2018-12-13

前期一篇文章:Modern GMT Series：Slice in 3D View (三維切片圖)中提到了在科研作圖中會經常遇到三維作圖的問題，而且GMT可以做三維圖且匯出的圖片質量非常高！但是也提到了當前GMT版本中的三維繪圖存在漏洞。主要兩個問題：（1）切片位置錯亂；（2）三維文字映象。就這兩個問題極大的限制了GMT三維作圖的應用。因為最近論文裡面要做三維立體圖(地理座標+笛卡爾座標+海底地形+剖面切片等)，找了找也沒找到更好的解決方案，乾脆自己動手豐衣足食——修改GMT的原始碼。經過一天半的努力，已經比較完美的修復了這兩個漏洞。目前官方版本中並沒有修復此bug。下面就介紹一下修復漏洞的過程和最終效果！

先上一張BUG修復後的效果圖鑑賞一下

三維地形圖+剖面切片+線+文字+方向標+ModernGMTlogo

注：由於我的論文正在審稿中，不便上傳較為複雜的繪圖結果，暫用一個相對簡單的圖來代替。

簡單例子

問題描述：繪製一個三維立方體並在每個面上標註文字以及設定每個面具有不同顏色且半透明；水平旋轉45度，然後x軸和y軸標註座標標籤。

GMT官方版本執行的結果

GMT6.0三維繪圖結果

繪圖程式碼

# Test 3D View of GMT: simple example
# Zhikui Guo, 2018-12-09, GEOMAR, Germany

function preset()
{
    figname=Test_3d
    angle_view=-45/25
    width_fig_x=10
    width_fig_y=10
    width_fig_z=7
    # 透明度
    alpha_profile=50
    # range
    xmin=0
    xmax=100
    ymin=0
    ymax=50
    zmin=-20
    zmax=0
    xc=`echo $xmin $xmax | awk '{print $1+($2-$1)/2}'`
    yc=`echo $ymin $ymax | awk '{print $1+($2-$1)/2}'`
    zc=`echo $zmin $zmax | awk '{print $1+($2-$1)/2}'`
    len_z=`echo $zmin $zmax | awk '{print ($2-$1)}'`
}
# 1. apply theme
# . styles.sh
# MonokaiTheme
# 常用程式碼標頭檔案:繪製logo
# . stdafx.sh
# 範圍和顏色等設定
preset

# plot 
gmt begin $figname pdf,png
    gmt basemap -JX$width_fig_x/$width_fig_y -JZ$width_fig_z -R$xmin/$xmax/$ymin/$ymax/$zmin/$zmax -Ba -Bza+l"Z(m)" -BwsenZ+gred  -pz$angle_view/$zmin
    echo "$xc $yc zmin plane" | gmt pstext -JZ -p -F+f20p,Helvetica-Bold,blue=thinner,white+jCM+a190 -Dj0c/0c
    gmt basemap -JX$width_fig_x/$width_fig_y -JZ$width_fig_z -R$xmin/$xmax/$ymin/$ymax/$zmin/$zmax -Bwsenz+ggray  -pz$angle_view/$zmax
    echo "$xc $yc zmax plane" | gmt pstext -JZ -p -F+f20p,Helvetica-Bold,blue=thinner,white+jCM+a190 -Dj0c/0c
    gmt basemap -JX$width_fig_y/$width_fig_z -JZ$width_fig_x -R$ymin/$ymax/$zmin/$zmax/$xmin/$xmax -Ba  
[email protected]$alpha_profile -Bx+l"Y(m)" --MAP_FRAME_PEN=1,[email protected] --MAP_ANNOT_OFFSET=-0.2   -px$angle_view/$xminn
    echo "$yc $zc ymin plane" | gmt pstext -JZ -p -F+f20p,Helvetica-Bold,blue=thinner,white+jCM+a45 -Dj0c/0c
    gmt basemap -JX$width_fig_x/$width_fig_z -JZ$width_fig_y -R$xmin/$xmax/$zmin/$zmax/$ymin/$ymax -Ba  
[email protected]$alpha_profile -Bx+l"X(m)" --MAP_FRAME_PEN=1,[email protected] --MAP_ANNOT_OFFSET=-0.2  -py$angle_view/$ymax
    echo "$xc $zc xminn plane" | gmt pstext -JZ -p -F+f20p,Helvetica-Bold,blue=thinner,white+jCM+a-45 -Dj0c/0c
    # 使用logo
    # add_logo 0 0.5 moderngmt -grid=false 
gmt end
# open $figname.pdf
# rm tmp* gmt.conf gmt.history

這段程式碼不涉及什麼資料，直接可以複製執行的。這裡我將主題和logo註釋了，這是我自己定義的，因為你的電腦上肯定是沒有相關的檔案和函式。直接複製上面的程式碼應該是可以直接出圖的。注意這是Mac系統下的指令碼，也可以在Linux系統下使用，但是windows系統下需要將變數引用從$xxx改為%xxx%。

都有什麼BUG？

從上面的第二個圖中可以看出，官方版本的三維繪圖結果存在一下幾個問題：

切片位置錯亂：需要向y軸正方向平移一個量
文字出現映象：x軸和y軸的標註文字都變成了映象且有一個翻轉角度
普通文字映象：除了座標軸的標註文字外，剖面切片上的文字的角度也有錯亂和映象現象

如何修復BUG？

簡單介紹以上三個BUG是如何修復的，具體見修復後的原始碼(我的資源庫中有介紹如何獲取原始碼和可用軟體)。

切片位置錯亂問題

GMT繪圖是基於PostScript的，三維繪圖的座標旋轉理論見GMT三維座標旋轉理論。經測試發現問題出現在gmt_plot.c原始檔裡面的gmt_plane_perspective函式

a = b = c = d = e = f = 0.0;
if (plane < 0)            /* Reset to original matrix */
    PSL_command (PSL, "PSL_GPP setmatrix\n");
else {    /* New perspective plane: compute all derivatives and use full matrix */
    if (plane >= GMT_ZW) level = gmt_z_to_zz (GMT, level);    /* First convert world z coordinate to projected z coordinate */
    switch (plane % 3) {
        case GMT_X:    /* Constant x, Convert y,z to x',y' */
            a = GMT->current.proj.z_project.sin_az;
            b = -GMT->current.proj.z_project.cos_az * GMT->current.proj.z_project.sin_el;
            c = 0.0;
            d = GMT->current.proj.z_project.cos_el;
            e = GMT->current.proj.z_project.x_off - level * GMT->current.proj.z_project.cos_az;
            f = GMT->current.proj.z_project.y_off - level * GMT->current.proj.z_project.sin_az * GMT->current.proj.z_project.sin_el;
            break;
        case GMT_Y:    /* Constant y. Convert x,z to x',y' */
            a = -GMT->current.proj.z_project.cos_az;
            b = -GMT->current.proj.z_project.sin_az * GMT->current.proj.z_project.sin_el;
            c = 0.0;
            d = GMT->current.proj.z_project.cos_el;
            e = GMT->current.proj.z_project.x_off + level * GMT->current.proj.z_project.sin_az;
            f = GMT->current.proj.z_project.y_off - level * GMT->current.proj.z_project.cos_az * GMT->current.proj.z_project.sin_el;
            break;
        case GMT_Z:    /* Constant z. Convert x,y to x',y' */
            a = -GMT->current.proj.z_project.cos_az;
            b = -GMT->current.proj.z_project.sin_az * GMT->current.proj.z_project.sin_el;
            c = GMT->current.proj.z_project.sin_az;
            d = -GMT->current.proj.z_project.cos_az * GMT->current.proj.z_project.sin_el;
            e = GMT->current.proj.z_project.x_off;
            f = GMT->current.proj.z_project.y_off + level * GMT->current.proj.z_project.cos_el;
            break;
    }
    // printf("ix: %f iy: %f\n",PSL->internal.x2ix,PSL->internal.y2iy);
    /* First restore the old matrix or save the old one when that was not done before */
    PSL_command (PSL, "%s [%g %g %g %g %g %g] concat\n",
        (GMT->current.proj.z_project.plane >= 0) ? "PSL_GPP setmatrix" : "/PSL_GPP matrix currentmatrix def",
        a, b, c, d, e * PSL->internal.x2ix, f * PSL->internal.y2iy);
}

正如GMT三維座標旋轉理論所提到的公式，gmt_plane_perspective函式的這段程式碼目的就是根據MGT繪圖命令中的-R, -J, -JZ, -p引數計算座標旋轉矩陣引數a,b,c,d,e,f，最後使用PSL_command將座標旋轉矩陣寫入ps檔案就能得到三維的效果。

GMT_Z的情況都是OK的，正如上圖所示三維座標軸主框架沒錯。但是 GMT_X和 GMT_Y的情況，也就是分別於X軸和Y軸垂直的切片的情況，出現了問題。經測試發現，問題在於e, f計算錯誤。三種情況的e, f值需要一致才不會發生未知平移錯亂。

解決方案

這裡有一個很容易的解決方案就是在case GMT_Z計算得到e,f值的時候，開啟一個臨時檔案將這兩個值儲存到臨時檔案裡面，然後如果執行到case GMT_X或者case GMT_Y的時候將這個臨時檔案裡面的e,f值讀入覆蓋掉程式計算的結果。這雖然不是一種完美的解決方案，但是可以實實在在的修復這個bug。以後有時間了仔細研究GMT內部是如何計算這些引數然後找到計算錯誤的位置，再做完美修改方案。修改原始碼後，重新編譯得到更新的gmt程式，然後再執行這個例子，得到的結果如下：

修復剖面位置問題後的結果

座標軸標註文字映象問題

上面一步修復了剖面位置問題，剖面座標軸標註文字依然存在映象翻轉的現象。PS的映象翻轉命令是-1 1 scale，如果相對某個文字進行映象翻轉，只需要在ps檔案裡面找到這個文字對應的程式碼，然後在前面加上這個命令即可實現映象翻轉。找到這個magic命令真是不容易，參見一個論壇。

依照這個邏輯，在原始碼中跟蹤找到繪製座標軸label的位置：gmt_plot.c檔案中的gmt_xy_axis函式，關鍵程式碼如下：

/* Finally do axis label */
if (A->label[0] && annotate && !gmt_M_axis_is_geo_strict (GMT, axis)) {
    unsigned int far_ = !below;
    char *this_label = (far_ && A->secondary_label[0]) ? A->secondary_label : A->label;    /* Get primary or secondary axis label */
    if (!MM_set) PSL_command (PSL, "/MM {%s%sM} def\n", neg ? "neg " : "", (axis != GMT_X) ? "exch " : "");
    form = gmt_setfont (GMT, &GMT->current.setting.font_label);
    PSL_command (PSL, "/PSL_LH ");
    PSL_deftextdim (PSL, "-h", GMT->current.setting.font_label.size, "M");
    PSL_command (PSL, "def\n");
    PSL_command (PSL, "/PSL_L_y PSL_A0_y PSL_A1_y mx %d add %sdef\n", PSL_IZ (PSL, GMT->current.setting.map_label_offset), (neg == horizontal) ? "PSL_LH add " : "");
    /* Move to new anchor point */
    PSL_command (PSL, "%d PSL_L_y MM\n", PSL_IZ (PSL, 0.5 * length));
    
    if (axis == GMT_Y && A->label_mode) {
        i = (below) ? PSL_MR : PSL_ML;
        PSL_plottext (PSL, 0.0, 0.0, -GMT->current.setting.font_label.size, this_label, 0.0, i, form);
    }
    else
        PSL_plottext (PSL, 0.0, 0.0, -GMT->current.setting.font_label.size, this_label, horizontal ? 0.0 : 90.0, PSL_BC, form);
    
}

程式碼段裡面呼叫PSL_plottext繪製座標軸label(也就是圖中的X(m)和Y(m))，為了方便管理和程式碼重用，我在gmt_plot.c裡面自定義了兩個函式: AxisLable_Flip_GMT_X_Y 和AxisTickLabel_Flip_GMT_X_Y分別對axis label和axis tick label進行操作。注意：tick label除了映象問題還存在一個180度的旋轉角度問題，這些都需要根據-p引數進行重新計算和判斷然後修復。

加入自定義函式後的程式碼段如下：

/* Move to new anchor point */
PSL_command (PSL, "%d PSL_L_y MM\n", PSL_IZ (PSL, 0.5 * length));

AxisLable_Flip_GMT_X_Y(GMT);// X labbel: 官方版本中的對於x和y剖面切片的文字出現了映象，所以這裡只對x和y剖面進行翻轉
if (axis == GMT_Y && A->label_mode) {
    i = (below) ? PSL_MR : PSL_ML;
    PSL_plottext (PSL, 0.0, 0.0, -GMT->current.setting.font_label.size, this_label, 0.0, i, form);
}
else
    PSL_plottext (PSL, 0.0, 0.0, -GMT->current.setting.font_label.size, this_label, horizontal ? 0.0 : 90.0, PSL_BC, form);

AxisLable_Flip_GMT_X_Y(GMT);//恢復翻轉

ticklabel 和 axis label都在同一個函式裡面，呼叫自定義函式AxisTickLabel_Flip_GMT_X_Y 修復之後的程式碼片段如下：

for (i = 0; i < nx1; i++) {
    if (gmtlib_annot_pos (GMT, val0, val1, T, &knots[i], &t_use)) continue;            /* Outside range */
    if (axis == GMT_Z && fabs (knots[i] - GMT->current.proj.z_level) < GMT_CONV8_LIMIT) continue;    /* Skip z annotation coinciding with z-level plane */
    if (GMT->current.setting.map_frame_type & GMT_IS_INSIDE && (fabs (knots[i] - val0) < GMT_CONV8_LIMIT || fabs (knots[i] - val1) < GMT_CONV8_LIMIT)) continue;    /* Skip annotation on edges when MAP_FRAME_TYPE = inside */
    if (!is_interval && plot_skip_second_annot (k, knots[i], knots_p, np, primary)) continue;    /* Secondary annotation skipped when coinciding with primary annotation */
    x = (*xyz_fwd) (GMT, t_use);    /* Convert to inches on the page */
    /* Move to new anchor point */
    PSL_command (PSL, "%d PSL_A%d_y MM\n", PSL_IZ (PSL, x), annot_pos);
    if (label_c && label_c[i] && label_c[i][0])
        strncpy (string, label_c[i], GMT_LEN256-1);
    else
        gmtlib_get_coordinate_label (GMT, string, &GMT->current.plot.calclock, format, T, knots[i]);    /* Get annotation string */
    double angle_text2=AxisTickLabel_Flip_GMT_X_Y(GMT,text_angle); //座標軸標註文字翻轉
    PSL_plottext (PSL, 0.0, 0.0, -font.size, string, angle_text2, justify, form);
    AxisTickLabel_Flip_GMT_X_Y(GMT,text_angle); //恢復：座標軸標註文字翻轉
}
if (!faro) PSL_command (PSL, "/PSL_A%d_y PSL_A%d_y PSL_AH%d add def\n", annot_pos, annot_pos, annot_pos);

這裡省略細節，感興趣的可以直接看我的原始碼

第二個BUG修復之後的結果

座標軸label和tick label映象問題修復後的效果

剖面文字的映象和旋轉角度問題

第二個BUG修復之後，座標軸上的文字正確歸為了。但是剖面上的文字還存在問題：文字映象以及某些情況下存在一個角度偏轉。找到繪製文字的程式碼位置：pstext.c檔案中的GMT_pstext函式，程式碼太長，這裡只貼關鍵程式碼附近的幾行：

        c_txt[m] = strdup (curr_txt);
        c_x[m] = plot_x;
        c_y[m] = plot_y;
        c_just[m] = T.block_justify;
        c_font[m] = T.font;
        m++;
    }
    else {
        PSL_plottext (PSL, plot_x, plot_y, T.font.size, curr_txt, T.paragraph_angle, T.block_justify, fmode);
    }
    if (Ctrl->A.active) T.paragraph_angle = save_angle;    /* Restore original angle */
}

就是在這裡呼叫PSL_plottext 函式進行文字繪製，而PSL_plottext 函式定義在postscriptlight.h裡面，函式實現在postscriptlight.c裡面。進入這兩個檔案裡面發現不太容易傳遞GMT_CTRL *GMT引數(GMT這個資料結構中包含了很多很多資訊)，為了相容性。重新自定義一個int PSL_plottext_mirror (int isMirrow,struct PSL_CTRL *PSL, .....)的函式，在這個函式中增加一個引數isMirrow。用PSL_plottext_mirror替換上面程式碼段中的PSL_plottext。與前面的幾個問題的修復方法類似，在pstext.c中新增函式int Text_Flip_GMT_X_Y(struct GMT_CTRL *GMT)判斷何時需要對文字進行映象翻轉，然後將這個引數傳遞給新定義的文字繪製函式PSL_plottext_mirror 去執行相應的操作。

c_txt[m] = strdup (curr_txt);
        c_x[m] = plot_x;
        c_y[m] = plot_y;
        c_just[m] = T.block_justify;
        c_font[m] = T.font;
        m++;
    }
    else {
        // PSL_command(PSL,"%%繪製文字主控函式\n");
        // Text_Flip_GMT_X_Y(GMT); //文字映象翻轉
        PSL_plottext_mirror (Text_Flip_GMT_X_Y(GMT),PSL, plot_x, plot_y, T.font.size, curr_txt, T.paragraph_angle, T.block_justify, fmode);
        // Text_Flip_GMT_X_Y(GMT); //恢復
    }
    if (Ctrl->A.active) T.paragraph_angle = save_angle;    /* Restore original angle */
}

使用自定義的文字繪製函式PSL_plottext_mirror 。

至此，上面三個問題已全部修復完畢，修復bug之後的GMT姑且稱之為ModernGMT(GMT私房菜)。執行結果如下。

ModernGMT執行結果

本人修復BUG之後的版本稱之為ModernGMT，這是一個持續更新的私房版本，不僅修復各種bug還會新增一些有用的地球物理重磁位場相關的程式。而且我每個星期都會同步官方更新，所以ModernGMT總是比官方版本新一些且功能多一些。但是由於我時間有限不能暫時還沒有把自己做的更新推送給官方團隊，以後會考慮！

ModernGMT三維繪圖結果

後記

經過上面的修復之後，目前繪製正常的三維圖應該是沒有什麼大的問題了。如果有人嘗試一些奇怪的想法，比如旋轉角度很畸形或者其他的引數不和常理，不保證不出問題。

依然存在的問題

當使用-JM投影方式的時候，如果資料範圍很大，比如100度X80度的跨度，由於投影計算過程中的誤差等原因，可能會出現切片不貼合的現象。如下圖所示：

大範圍三維繪圖切片長度計算誤差

解決這個問題的技巧就是直接在GMT繪圖程式碼裡面調整一下切片的長度即可，比如增加一個小的倍數，稍微調整然後看著沒問題就行。以後有空了會從程式碼裡面直接解決

    # # profile along lat
    width_fig_y0=`echo $width_fig_y | awk '{print $1+0.3}'`
    gmt psbasemap -R$range_LatZLon -JX$width_fig_y/$width_fig_z -JZ$width_fig_x -BS -Ba --MAP_ANNOT_OFFSET=$offset_ticklabel_lon  -px$angle_view/$pos_profile_lon

小範圍是沒有問題的，比如下面這張圖

西南印度洋脊龍旗熱液區附近的ETOPO1水深資料

兩個剖面上的資料是隨便用程式生成的高斯分佈資料，不代表任何意義，只是為了說明在剖面切片中繪製grdimage的問題。

GMT更新程式獲取方法

由於時間和精力有限，暫不為伸手黨直接開放原始碼和程式。原始碼以及編譯之後的安裝檔案我會上傳到雲端，請到我的資源庫中檢視程式碼和程式獲取方法。

GMT三維繪圖有BUG? 修復它！

簡單例子

GMT官方版本執行的結果

繪圖程式碼

都有什麼BUG？

如何修復BUG？

切片位置錯亂問題

解決方案

座標軸標註文字映象問題

剖面文字的映象和旋轉角度問題

ModernGMT執行結果

後記

依然存在的問題

GMT更新程式獲取方法

GMT三維繪圖有BUG? 修復它！

45.Qt openGL實現三維繪圖

matlab的三維繪圖和四維繪圖

Python三維繪圖--Matplotlib

MATLAB中surf與scatter3不能同時三維繪圖解決辦法

Qt學習：三維繪圖之OpenGL和Qt的結合

Matplotlib Toolkits：三維繪圖工具包matplotlib.mplot3d

MFC+OpenGL三維繪圖(二)——開啟一個STL檔案並顯示

要怎麽在CAD編輯器中繪制三維線段？有什麽技巧？

CAD繪圖設定中怎麼進行設定三維建模?怎麼進行檢視有效模型？

C++ windows已在xxx.exe中觸發一個斷點,其原因可能是堆被損壞,這說明xx.exe中或它所載入的任何DLL中有bug。

CAD三維設計繪圖軟體AutoCAD 2017 for Mac破解教程

VC中“其原因可能是堆被損壞，這也說明*.exe中或它所載入的任何DLL中有bug”的分析

2018華為杯E題的一些思考和三維模型圖繪圖技巧總結

(視訊)《快速建立網站》 4.1 為啥造軟體不同於造汽車，為啥是軟體就一定會有Bug －構建開發運維一體化（DevOps)

其原因可能是堆被損壞，這說明**.exe中或它載入的任何DLL中有Bug

機房三維（3D）監控系統和機房視覺化動力環境監控系統，兩者有什麼功能特點？

暫存，本人博客有bug，正在全力修復。

”其原因可能是堆被損壞，這也說明exe中或它所載入的任何DLL 中有bug“的解決辦法

【Bug修復】toolbar頂端與狀態列有間隙

GMT三維繪圖有BUG? 修復它！

簡單例子

GMT官方版本執行的結果

繪圖程式碼

都有什麼BUG？

如何修復BUG？

切片位置錯亂問題

解決方案

座標軸標註文字映象問題

剖面文字的映象和旋轉角度問題

ModernGMT執行結果

後記

依然存在的問題

GMT更新程式獲取方法

相關推薦