球坐标中等矩网格数据转换的递推算法

球坐标中等矩网格数据转换的递推算法 球坐标中等矩网格数据转换的递推算法 陈维海 唐泽圣 ( )清华大学计算机系 北京 100084 摘要 本文为解决气象领域可视化中所遇到的网格数据转换问题, 提出了等矩网格数据转换 的递推算法, 并分析了算法的复杂性, 本文最后给出了算法在交互式可视化系统M e teoV is 中的应用结果。 关 键 词 递推算法, 算法复杂性, 可视化。 1 引 言 在科学计算可视化系统中, 网格数据及其赖以存在的坐标系, 都随着应用领域的变化 而变化。而网格数据在不同坐标系中的转换, 亦是可视化系统经常遇到的问题。在气象数 1 据的可视化系统中, 原始数据通常在等间矩的经纬网格及不等间矩的高度网 M e teoV is 格中测定与计算。 将数据映射到三维空间时, 通常采用两种方法: ( ) 1局域坐标系: 将地球表面近似为平面, 将经线、纬线看作两个相互垂直的坐标方 向, 与高度一起, 形成三维直角坐标系。这种方法简便易行, 在中纬度地区亦能较精确地反 2- 4 映数据内涵, 因而被许多系统采用。 然而, 当原始数据范围增大, 所处纬度偏高时, 球 面效应引起的误差已不能由于平面近似而忽略不计, 从而产生了第二种映射方法。 () 2球坐标系: 将地球近似为圆球, 在球坐标中显示观测及计算数据。 此方法能够最 真实地反映气象数据本身的含义, 因而深受气象工作者的欢迎。 然而, 在球坐标系的显示 中, 由于从球坐标到直角坐标的数据转换引入了大量超越函数的计算, 使得系统整体性能 大大降低。 本文意在利用气象数据网格分布的特点, 提出一套转换算法的递推公式, 以解 决上面所述的问题。 2 记号约定 为叙述及书写方便, 在本文中特作下述约定: 记 C = co s? Α, S = sin ? Α, C = co s? Β, S = sin ? Β,Α Α Β Β m m m 11 12 13 C 矩阵M = ; 表示超越运算 sin 或 co s; ? 表示乘法运算; ? 表示加法 m m m 21 22 23 m m m 31 32 33 运算。 3 球坐标系中的坐标转换 () 如图 1 所示, 已知球坐标 Α, Β, r, 则其对应的直角坐标公式 1 为: x = rco sΒco sΑ y = rco sΒsin Α z = r sin Β 根据等矩网格数据的 特 点: ? Α, ? Β, ? r 恒 定, 可分别得到三组公式。 () 3. 1 Α 改变时的递推公式 公式 2 () Α, Β, r, Α递增值为 ? Α,设初始坐标为 0 0 0 则: () ()1根据公式 1, 计算初始坐标 x , y , z 0 0 0 x = rco sΒco sΑ0 0 0 y = rco sΒsin Α0 0 0 图 1 z = r sin Β0 0 计算 C Α= co s ? Α, S Α= sin ? Α () () 2对于每一步, 设 x ′, y ′, z ′为递推计算后的新坐标, 则: x ′= x r C - y r S Α Α y ′= x S + y C r r Α Α z ′= z 推导: 将 Α= Α′+ ? Α代入公式 1 得 x ′= rco sΒco sΑ′ ()= rco sΒco sΑco s? Α- sin Αsin ? Α = rco sΒco sΑco s? Α- rco sΒsin Αsin ? Α = x co s? Α- y sin ? Α y ′= rco sΑsin Α′ ( )= rco sΒsin Αco s? Α+ co sΑsin ? Α = rco sΒsin Αco s? Α+ rco sΒco sΑsin ? Α = y co s? Α+ x sin ? Α z ′= r sin Β = z (设初始坐) 3. 2 Β 改变时的递推公式 公式 3 () () 标为 Α, Β, r, Β 递增值为 ? Β: 1根据0 0 0 ()公式 1 计算初始坐标 x , y , z 0 0 0 x = rco sΒco sΑ0 0 0 y = rco sΒsin Α0 0 0 z = r sin Β0 0 计算 C = co s? Β, S = sin ? Β, 令 K = S co sΑ, K = S sin Α, K = S ƒco sΑΒ Β 1 Β2 Β3 Β () 2递推公式为: x ′= C r x - K r z Β1 y ′= C y - K z rr Β2 z ′= C r z - K r xΒ3 推导: ) (x ′= rco s Β + ? Βco sΑ = rco sΒco s? Βco sΑ- r sin Βsin ? Βco sΑ = x co s? Β - z co sΑsin ? Β 因为 co sΑ为常数, 所以 x ′= C - K zΒ 1 ) (y ′= rco s Β + ? Βsin Α = rco sΒco s? Βsin Α- r sin Βsin ? Βsin Α = x sin ? Β - z sin Αsin ? Β 因为 co sΑ为常数, 所以 y ′= C r y - K z Β2 )(z ′= r sin Β + ? Β = r sin Βco s? Β + rco sΒsin ? Β ( ) ( )= z co s? Β + rco sΒco sΑsin ? Βƒco sΑ = C + K xΒ 3 (设初始坐标为 ) 3. 3 r 改变时的递推公式 公式 4 () () Α, Β, r, r 的递增值为 ? r 则: 1应用公式 0 0 0 () 1 求出 x , y , z , 并令0 0 0 ? x = ? rco sΒco sΑ, ? y = ? rco sΒsin Α, ? z = ? r sin Β () 2递推公式: x ′= x + ? x , y ′= y + ? y , z ′= z + ? z推导过程略。 3. 4 算法分析 算法分析均以m × n 水平网格为例, 其中m 为东西向的网格数目, n 为南北向网格数 目。 3. 4. 1 应用公式 1 C C 由于每转换一组坐标需进行 5 ? + 3 , 所以m × n 网格需运算为 3m n ? + 5m n , () 即乘法与超越运算次数均为 O m n 级。 3. 4. 2 应用公式 2 及应用公式 3 C 应用公式 2, 需进行初始计算 5 + 5 ? , 而每一步递推, 需计算 4 ? + 2 ? 。 C 应用公式 3, 需进行初始计算 7 + 8 ? , 而每一步递推, 需计算 6 ? + 3 ? 。 由于 ? Α, ? Β 均为常数, 因而可将 m × n 的转换算法优化如下: () 1计算 C , S ΑΒ 公式 3 初始计算 () 2 fo r j = 1 to n do () 3() ()输出 x ′, y ′, z ′坐标 此时已完成公式 2 的初始计算 () 4 fo r i = 2 to m do () 5应用公式 2 () 6() 输出 x ′, y ′, z ′坐标 () 7end do () 8if j < n th en 应用公式 3 ()9 end do ()10 C C () () (() )()其中 5— 8步共需 m - 1〔4 ? + 2 ? 〕次计算, 1、2步共需 2 + 〔7 + 8 ? 〕次计算。因此完成 m × n 次转换共需计算为: C )() (9 + 8 ? + n m - 1〔4 ? + 2 ? 〕+ n - 1〔6 ? + 3 ? 〕 C = m n〔4 ? + 2 ? 〕+ n〔6 ? + 3 ? - 4 ? - 2 ? 〕+ 9 + 8 ? - 6 ? - 3 ? C = m n〔4 ? + 2 ? 〕+ n〔2 ? + ? 〕+ 9 + 2 ? - 3 ? 与 3. 4. 1 的结论相比, 递推算法将原公式的 3m n 次超越运算与 5m n 次乘法运算减少 为约 4m n 次乘法及 2m n 次加法, 由于超越运算的计算时间大大多于乘法运算, 因而应用 递推公式能够大大提高转换效率。 4 球坐标中的矢量转换 ?() 如图 2 所示, 处于 Α, Β, r的矢量 v= u v w , 可通过下列矩阵变换到直角坐标系 ?中的矢量 v ′= x y z : ? ? v ′= vM r 其中公式 5:0 0 () () () ()M = R o t90r R o t90r R o t- Βr R o tΑz y y z - sin Α co sΑ 0 = - co sΑsin Β - sin Αsin Β co sΒ - co sΑco sΒ sin Αco sΒ sin Β 根据网格的等矩特性, 亦可得到以下两组矩阵 递推公式。 ()4. 1 当 Α 改变时的递推公式 公式 6 () 设 矢 量 的 初 始 坐 标 为 Α, Β, r, Α 递 增 值 为0 0 0 ? Α, 则: 图 2 () 1应用公式 5 计算M 0 Αco sΑ0- sin 0 0 - co sΑsin Β- sin Αsin Βco sΒ0 0 0 0 0 M =0 co sΑco sΒsin Αco sΒsin Β0 0 0 0 0 计算 C Α= co s ? Α, S Α= sin ? Α () 2对于每一步, 设M ′为递推计算后的新矩阵, 则 () () m C - m S m C + m S 011Α 12Α12Α 11Α M ′= () () m C - m S m C + m S m 21Α 22Α22Α 21Α23 () () m C - m S m C + m S m 31Α 32Α32Α 31Α33 推导略。 (设矢量的初始坐) 4. 2 当 Β 改变时的递推公式 公式 7 () () 标为 Α, Β, r, Β 递增值为 ? Β, 则: 1应用公式 5 0 0 0 计算M 0 Αco sΑ0- sin 0 0 - co sΑsin Β- sin Αsin Βco sΒ0 0 0 0 0 M =0 co sΑco sΒsin Αco sΒsin Β0 0 0 0 0 计算 C = co s? Β, S = sin ? ΒΒ Β () 2对于每一步, 设M ′为递推计算后的新矩阵, 则: m 11 m 12 0 M ′= () () () m C - m S m C - m S m C - m S 21Β 31Β22Β 32Β23Β 33Β () () ()m C + m S m C + m S m C + m S 21Β 31Β22Β 32Β23Β 33Β 推导略。 4. 3 算法分析 与 3. 4 所使用的网格一致, 此处也使用 m × n 水平网格矢量。 4. 3. 1 应用公式 5 C C 由于每转换一组坐标需进行 4 ? + 4 , 所以m ×n 网格需运算 4m n ? + 4m n , 即乘 () 法与超越运算次数均为 O m n 级。 4. 3. 2 应用公式 6 及应用公式 7 C 应用公式 6, 需进行初始计算 6 + 4 ? , 而每一步递推, 需计算 12 ? + 6 ? 。 C 应用公式 7, 需进行初始计算 6 + 4 ? , 而每一步递推, 需计算 12 ? + 6 ? 。 由于 ? Α, ? Β 均为常数, 因而可将 m × n 的转换算法优化如下: 计算 C , S ΑΒ () 1公式 7 初始计算 () 2 fo r j = 1 to n do () 3()输出M ′此时已完成公式 6 的初始计算 () 4fo r i = 2 to m do () 5应用公式 6 () 6输出M ′ () 7end do () 8if j < n th en 应用公式 7 ()9 end do ()10 C C () () (() )()其中 5— 8步共需 m - 1〔12 ? + 6 ? 〕次计算, 1、2步共需 2 + 〔6 + 4 ? 〕次计 算。因此完成 m × n 次转换共需计算为: C ()()8 + 4 ? + n m - 1〔12 ? + 6 ? 〕+ n - 1〔12 ? + 6 ? 〕 C ()= m n - 1〔12 ? + 6 ? 〕+ 8 + 4 ? C = m n〔12 ? + 6 ? 〕+ 8 - 8 ? - 6 ? 与 4. 3. 1 中的结论相比, 由于在计算机中超越函数的计算时间远远大于乘法运算的 时间, 因而此递推算法亦能提高运算速度。 5 实现与结论 本文所述的递推公式, 均已在气象数据可视化系统M e teoV is 中得以实现。图 3 即为 水平风矢量在局域坐标系中的显示效果, 图 4 为同一数据在球坐标系中的显示效果。图中 箭头所指的方向为风矢量的方向, 箭头的长度为风矢量的模长。由此可以看出两者的差距 以及球坐标显示的必要性。 图 3 图 4 由于是一个高度交互的气象可视化系统, 因此, 使用递推公式处理球坐标 M e teoV is 下的显示模型, 大大加快了显示速度, 提高了系统的交互性能。 参 考 文 献 1 ] , . . ’94, C h en W H T ang Z SA h igh ly in te rac t ive m e teo ro lo g ica l v isua liza t io n sy stem M e teoV isCA DDM A ug 1994, 117, 122. 2 ] , . . 22: 1988,Sch iavo ne J A J u le sz BA pp lica t io n o f com p u te r g rap h ic s to th e v isua liza t io n o f m e teo ro lo g ica l da ta 327, 334. 3 ] , . . , Ge lbe rg L M S tep h en so n T PSup e rcom p u t ing g rap h ic s in ea r th and p lane ta ry sc ienceIE E E C G&A J u ly 1987, 23, 26. 4 ] . . , 1992, 10, 12.T re in ish L ACo r re la t ive da ta ana ly sis in th e ea r th sc ienceIE E E C G&A M ay A N ITERA T IVE AL GO R ITHM FO R CO NVERS IO N O F - ISO GR ID D A TA IN SPHER ICAL COO RD INA TES C h en W e ih a i an d T an g Ze sh en g (). . . . , , 100084D ep tof C om pS ci& T echT sing h u a U n iv ersity B eij ing , A bstra c t In o rde r to co n ve r t th e g r id da ta in v isu a liza t io n o f m e teo ro lo g ica l da taan 2. 2 ite ra t ive a lgo r ithm fo r co n ve r sio n o f iso g r id da ta is p re sen tedT h e a lgo r ithm com p lex i ty is an a lyzed an d th e app lica t io n re su lt in a n ew in te rac t ive m e teo ro lo g ica l v isu a liza t io n 2.sy stem M e teoV is is show n , , . Key word s ite ra t ive a lgo r ithm a lgo r ithm com p lex ityv isu a liza t io n 第五届国际计算机辅助 设计 领导形象设计圆作业设计ao工艺污水处理厂设计附属工程施工组织设计清扫机器人结构设计 及图形学学术会议 ) ( 第一轮通知ƒ97CAD Graph ic s’ 会议日期: 1997 年 11 月底—12 月初 会议地点: 深圳 主办单位: 中国计算机学会, 中国科学院计算技术研究所 开放实验室 CA D 征文截止日期: 1997 年 4 月 30 日 () 征文范围 但不限于: ? 计算机辅助设计 ? 计算机图形学 ? 工程数据库 ? 科学计算可视化 ? CA D ƒCAM 技术 ? 真实感图形绘制技术 ? ? 虚拟现实 C IM S ? 计算机支持的协同设计 ? 自然景物模拟 ? 计算机辅助几何设计 ? 多媒体设计 ? 测试、论断与可测性设计 ? 计算机动画技术 ? 容错技术 ? 人机交互技术 论文交寄地址:100080 北京 2704 信箱 开放实验室 李 华收 CA D 电话 + 8610 62565533—5658 论文要求: () 11 论文须用英文书写; 字数 包括插图不超过 6000 字 ( )( )21 论文除正文外, 还必须包括论文摘要 不超过 200 字、关键字 最多 8 个字、作者姓名、所在单 位、通信地址、电话号码及下面文字“, if th e p ap e r is accep tedo ne o f th e au tho r s w ill p re sen t th e p ap e r a t ”和个人签名; 如有可能, 请附上 2地址及 号码; 通信地址包括 ƒ97 CA D G rap h ic s’co nfe renceE m a il FA X 中、英文。 31 论文必须是没有发表过的; 41 论文一式四份。 欢迎踊跃投稿。