Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
Browse files

Added utility files to rasterize the mesh

  • Loading branch information...
commit a40925a1335756538201c04812b97e4e6134d12b 1 parent 29535f3
Chris Marsh authored
128 utils/rasterive_mesh/SaveAsciiRaster.m
View
@@ -0,0 +1,128 @@
+function SaveAsciiRaster(varname, header,fname);
+ % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
+ % %
+ % Produced by Giuliano Langella %
+ % e-mail:gyuliano@libero.it %
+ % March 2008 %
+ % %
+ % Last Updated: 01 June, 2009 %
+ % %
+ % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
+
+
+%§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§
+%
+% ----- SYNTAX -----
+% SaveAsciiRaster(varname, header);
+% SaveAsciiRaster(varname); % in case varname is an xyz matrix
+%
+%
+% ----- DESCRIPTION -----
+% This function saves a spatial matrix into an Arc-Info ascii raster. Two file extension '.asc'
+% or '.txt' are supported.
+% FIRST CASE
+% USE: SaveAsciiRaster(varname, header);
+% It requires two inputs: (1) the z-values to be exported ('varname'
+% variable), and (2) the 'header' vector with the spatial information of
+% the grid. 'varname' can be a 1-D vector or a 2-D spatial grid.
+% SECOND CASE
+% USE: SaveAsciiRaster(varname);
+% If an xyz matrix (with [x_coord,y_coord,z_values]) is given as
+% 'varname', no 'header' has to be defined, since the function will
+% extract all the required header information from the xyz table. The
+% first row contains the x_coord, y_coord and z_value of the most
+% north-western cell; the last row refers to the most south-eastern
+% pixel. Elements in xyz are sorted column-by-column from the
+% geographical grid
+% (geographical_grid=[1st_col,2nd_col,3rd_col,...,last_col];
+% xyz=[1st_col;2nd_col;3rd_col;...;last_col]).
+% The xy coordinates have to refer to the center of the cells.
+%
+%
+% ----- INPUT -----
+% +'varname'(mandatory) : the MatLab matrix to be exported in ascii raster.
+% It can be passed:
+% (1) a 2-D z-values matrix of size equal to the spatial grid extent;
+% (2) a 1-D vector with z-values sorted as [1st_col;2nd_col;...];
+% (3) an xyz table with characteristics described in "SECOND CASE".
+%
+%
+% +'header'(facultative): the geospatial reference matrix.
+% If varname is an xyz matrix then header have not to be given.
+% The header matrix as the .hdr file in Arc-Info Binary Raster.
+% The 'header' MatLab variable is created when importing an ascii
+% raster with function ImportAsciiRaster (author: Giuliano Langella).
+%
+%
+% ----- EXAMPLE1 -----
+% [Z h] = ImportAsciiRaster(NaN, 'h'); % import DEM and Aspect
+% filt_gau05 = fspecial('gaussian', 3, 0.5); % create a low-pass spatial filter
+% Z_g5 = imfilter(Z, filt_gau05, 'symmetric'); % apply the filter
+% SaveAsciiRaster(Z_g5, h); % save the filtered DEM
+% ----- EXAMPLE2 -----
+% SaveAsciiRaster(xyz); % save the xyz matrix as an ascii grid
+%§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§
+
+
+
+% Save path
+% [FileName, PathName] = uiputfile({'*.asc','Arc-Info ASCII Raster (*.asc)'; ...
+% '*.txt','Text ASCII Raster (*.txt)'; ...
+% '*.*', 'All Files (*.*)'}, 'Save Ascii Grid', pwd);
+
+% OPEN file
+fid = fopen(fname,'w');
+
+% if I loaded an xyz matrix as 'varname'
+if nargin == 1
+ cellsize = abs(varname(1,2)-varname(2,2));
+ hor = max(varname(:,1)) - min(varname(:,1)) + cellsize;
+ ver = max(varname(:,2)) - min(varname(:,2)) + cellsize;
+ ncols = ceil(hor/cellsize);
+ nrows = ceil(ver/cellsize);
+ raster = zeros(nrows,ncols);
+ %create the header variable [default no data value is -9999]
+ header = [ncols; nrows; min(varname(:,1))-0.5*cellsize; min(varname(:,2))-0.5*cellsize; cellsize; -9999];
+ Zvar = reshape(varname(:,3),header(2),header(1));
+ varname=[];
+ varname=Zvar;
+end
+
+% WRITE HEADER
+fprintf(fid,'%s','ncols '); %1
+fprintf(fid,'%12.0f\n', header(1,1));
+fprintf(fid,'%s','nrows '); %2
+fprintf(fid,'%12.0f\n', header(2,1));
+fprintf(fid,'%s','xllcorner '); %3
+fprintf(fid,'%f\n', header(3,1));
+fprintf(fid,'%s','yllcorner '); %4
+fprintf(fid,'%f\n', header(4,1));
+fprintf(fid,'%s','cellsize '); %5
+fprintf(fid,'%f\n', header(5,1));
+fprintf(fid,'%s','NODATA_value '); %6
+fprintf(fid,'%f\n', header(6,1));
+
+% WRITE MATRIX
+%substitute to NaN the NODATA_value written in header
+varname(find(isnan(varname))) = header(6,1);
+%start loop
+ncols = header(1,1);
+nrows = header(2,1);
+handle = waitbar(0, mfilename);
+for CurrRow = 1:nrows;
+ waitbar(CurrRow/nrows);
+ % if varname is a vector instead of a 2-D array
+ if size(varname,2) == 1;
+ fprintf(fid,'% f ',varname( ((CurrRow-1)*ncols + 1) : (CurrRow*ncols) )' );
+ fprintf(fid,'%s\n', ' ');
+ % if varname is a 2-D array
+ else
+ fprintf(fid,'%f ',varname(CurrRow,:));
+ fprintf(fid,'%s\n', ' ');
+ end
+end
+fclose(fid);
+fclose('all');
+
+waitbar(CurrRow/nrows, handle, 'Done!');
+close(handle)
84 utils/rasterive_mesh/gridtrimesh.m
View
@@ -0,0 +1,84 @@
+function Z = gridtrimesh(F,V,X,Y,faceValue)
+%GRIDTRIMESH Fitting a square grid onto a triangular mesh.
+%
+% Z = GRIDTRIMESH(F,V,X,Y) fits a surface of the form Z = F(X,Y) to the
+% triangular mesh defined by F and V.
+%
+% The triangles of the triangular mesh are defined in the m-by-3 face
+% matrix F. Each row of F defines a single triangular face by indexing
+% into the n-by-3 matrix V that contains X, Y and Z coordinates of the
+% vertices.
+%
+% It is assumed that X and Y are produced by [X,Y] = meshgrid(x,y),
+% where x and y are monotonically increasing vectors.
+%
+% The function Z = MXGRIDTRIMESH(F,V,X,Y) gives exactly the same result
+% but is coded in C and compiled using MEX resulting in slightly faster
+% execution.
+%
+% Example: (assume the file example.mat contains F and V)
+% load example
+% nx = 10; ny = 10;
+% x = linspace(min(V(:,1)),max(V(:,1)),nx);
+% y = linspace(min(V(:,2)),max(V(:,2)),ny);
+% [X,Y] = meshgrid(x,y);
+% Z = gridtrimesh(F,V,X,Y);
+% surf(X,Y,Z)
+%
+%
+% Author: Willie Brink [ w.brink@shu.ac.uk ]
+% April 2007
+
+
+
+m = size(F,1);
+
+% size of grid
+nx = size(X,2);
+ny = size(X,1);
+
+% initialise output
+Z = NaN(size(X));
+
+% consider every triangle projected to the x-y plane and determine whether gridpoints lie inside
+for i = 1:m,
+ v1x = V(F(i,1),1); v1y = V(F(i,1),2); v1z = V(F(i,1),3);
+ v2x = V(F(i,2),1); v2y = V(F(i,2),2); v2z = V(F(i,2),3);
+ v3x = V(F(i,3),1); v3y = V(F(i,3),2); v3z = V(F(i,3),3);
+ % we'll use the projected triangle's bounding box: of the form (minx,maxx) x (miny,maxy)
+ minx = min([v1x v2x v3x]);
+ maxx = max([v1x v2x v3x]);
+ % find smallest x-grid value > minx, and largest x-grid value < maxx
+ east = NaN; west = NaN;
+ j = 1; while (j <= nx && X(1,j) < minx), j = j + 1; end; if (j <= nx), west = j; end
+ j = nx; while (j >= 1 && X(1,j) > maxx), j = j - 1; end; if (j >= 1), east = j; end
+ % if there are gridpoints strictly inside bounds (minx,maxx), continue
+ if ~isnan(east) && ~isnan(west) && east - west >= 0,
+ miny = min([v1y v2y v3y]);
+ maxy = max([v1y v2y v3y]);
+ % find smallest y-grid value > miny, and largest y-grid value < maxy
+ north = NaN; south = NaN;
+ j = 1; while (j <= ny && Y(j,1) < miny), j = j + 1; end; if (j <= ny), north = j; end
+ j = ny; while (j >= 1 && Y(j,1) > maxy), j = j - 1; end; if (j >= 0), south = j; end
+ % if, further, there are gridpoints strictly inside bounds (miny,maxy), continue
+ if ~isnan(north) && ~isnan(south) && south - north >= 0,
+ % we now know that there might be gridpoints bounded by (west,east) x (north,south)
+ % that lie inside the current triangle, so we'll test each of them
+ for j = west:east,
+ for k = north:south,
+ % calculate barycentric coordinates of gridpoint w.r.t. current (projected) triangle
+ A = [v1x v2x v3x; v1y v2y v3y; 1 1 1];
+ if rcond(A) > eps, w = A\[X(k,j); Y(k,j); 1]; else w = [1; 1; 1]/3; end
+ if min(w) > 0,
+ % use barycentric coordinates to calculate z-value
+ z = w(1)*v1z + w(2)*v2z + w(3)*v3z;
+ if isnan(Z(k,j)) || z > Z(k,j)
+ Z(k,j) = faceValue(i);
+% Z(k,j) = z;
+ end
+ end
+ end
+ end
+ end
+ end
+end
231 utils/rasterive_mesh/minboundrect.m
View
@@ -0,0 +1,231 @@
+% Copyright (c) 2009, John D'Errico
+% All rights reserved.
+%
+% Redistribution and use in source and binary forms, with or without
+% modification, are permitted provided that the following conditions are
+% met:
+%
+% * Redistributions of source code must retain the above copyright
+% notice, this list of conditions and the following disclaimer.
+% * Redistributions in binary form must reproduce the above copyright
+% notice, this list of conditions and the following disclaimer in
+% the documentation and/or other materials provided with the distribution
+%
+% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+% POSSIBILITY OF SUCH DAMAGE.
+% from http://www.mathworks.com/matlabcentral/fileexchange/13624-minimal-bounding-rectangle
+
+
+function [rectx,recty,area,perimeter] = minboundrect(x,y,metric)
+% minboundrect: Compute the minimal bounding rectangle of points in the plane
+% usage: [rectx,recty,area,perimeter] = minboundrect(x,y,metric)
+%
+% arguments: (input)
+% x,y - vectors of points, describing points in the plane as
+% (x,y) pairs. x and y must be the same lengths.
+%
+% metric - (OPTIONAL) - single letter character flag which
+% denotes the use of minimal area or perimeter as the
+% metric to be minimized. metric may be either 'a' or 'p',
+% capitalization is ignored. Any other contraction of 'area'
+% or 'perimeter' is also accepted.
+%
+% DEFAULT: 'a' ('area')
+%
+% arguments: (output)
+% rectx,recty - 5x1 vectors of points that define the minimal
+% bounding rectangle.
+%
+% area - (scalar) area of the minimal rect itself.
+%
+% perimeter - (scalar) perimeter of the minimal rect as found
+%
+%
+% Note: For those individuals who would prefer the rect with minimum
+% perimeter or area, careful testing convinces me that the minimum area
+% rect was generally also the minimum perimeter rect on most problems
+% (with one class of exceptions). This same testing appeared to verify my
+% assumption that the minimum area rect must always contain at least
+% one edge of the convex hull. The exception I refer to above is for
+% problems when the convex hull is composed of only a few points,
+% most likely exactly 3. Here one may see differences between the
+% two metrics. My thanks to Roger Stafford for pointing out this
+% class of counter-examples.
+%
+% Thanks are also due to Roger for pointing out a proof that the
+% bounding rect must always contain an edge of the convex hull, in
+% both the minimal perimeter and area cases.
+%
+%
+% Example usage:
+% x = rand(50000,1);
+% y = rand(50000,1);
+% tic,[rx,ry,area] = minboundrect(x,y);toc
+%
+% Elapsed time is 0.105754 seconds.
+%
+% [rx,ry]
+% ans =
+% 0.99994 -4.2515e-06
+% 0.99998 0.99999
+% 2.6441e-05 1
+% -5.1673e-06 2.7356e-05
+% 0.99994 -4.2515e-06
+%
+% area
+% area =
+% 0.99994
+%
+%
+% See also: minboundcircle, minboundtri, minboundsphere
+%
+%
+% Author: John D'Errico
+% E-mail: woodchips@rochester.rr.com
+% Release: 3.0
+% Release date: 3/7/07
+
+% default for metric
+if (nargin<3) || isempty(metric)
+ metric = 'a';
+elseif ~ischar(metric)
+ error 'metric must be a character flag if it is supplied.'
+else
+ % check for 'a' or 'p'
+ metric = lower(metric(:)');
+ ind = strmatch(metric,{'area','perimeter'});
+ if isempty(ind)
+ error 'metric does not match either ''area'' or ''perimeter'''
+ end
+
+ % just keep the first letter.
+ metric = metric(1);
+end
+
+% preprocess data
+x=x(:);
+y=y(:);
+
+% not many error checks to worry about
+n = length(x);
+if n~=length(y)
+ error 'x and y must be the same sizes'
+end
+
+% start out with the convex hull of the points to
+% reduce the problem dramatically. Note that any
+% points in the interior of the convex hull are
+% never needed, so we drop them.
+if n>3
+ edges = convhull(x,y); % 'Pp' will silence the warnings
+
+ % exclude those points inside the hull as not relevant
+ % also sorts the points into their convex hull as a
+ % closed polygon
+
+ x = x(edges);
+ y = y(edges);
+
+ % probably fewer points now, unless the points are fully convex
+ nedges = length(x) - 1;
+elseif n>1
+ % n must be 2 or 3
+ nedges = n;
+ x(end+1) = x(1);
+ y(end+1) = y(1);
+else
+ % n must be 0 or 1
+ nedges = n;
+end
+
+% now we must find the bounding rectangle of those
+% that remain.
+
+% special case small numbers of points. If we trip any
+% of these cases, then we are done, so return.
+switch nedges
+ case 0
+ % empty begets empty
+ rectx = [];
+ recty = [];
+ area = [];
+ perimeter = [];
+ return
+ case 1
+ % with one point, the rect is simple.
+ rectx = repmat(x,1,5);
+ recty = repmat(y,1,5);
+ area = 0;
+ perimeter = 0;
+ return
+ case 2
+ % only two points. also simple.
+ rectx = x([1 2 2 1 1]);
+ recty = y([1 2 2 1 1]);
+ area = 0;
+ perimeter = 2*sqrt(diff(x).^2 + diff(y).^2);
+ return
+end
+% 3 or more points.
+
+% will need a 2x2 rotation matrix through an angle theta
+Rmat = @(theta) [cos(theta) sin(theta);-sin(theta) cos(theta)];
+
+% get the angle of each edge of the hull polygon.
+ind = 1:(length(x)-1);
+edgeangles = atan2(y(ind+1) - y(ind),x(ind+1) - x(ind));
+% move the angle into the first quadrant.
+edgeangles = unique(mod(edgeangles,pi/2));
+
+% now just check each edge of the hull
+nang = length(edgeangles);
+area = inf;
+perimeter = inf;
+met = inf;
+xy = [x,y];
+for i = 1:nang
+ % rotate the data through -theta
+ rot = Rmat(-edgeangles(i));
+ xyr = xy*rot;
+ xymin = min(xyr,[],1);
+ xymax = max(xyr,[],1);
+
+ % The area is simple, as is the perimeter
+ A_i = prod(xymax - xymin);
+ P_i = 2*sum(xymax-xymin);
+
+ if metric=='a'
+ M_i = A_i;
+ else
+ M_i = P_i;
+ end
+
+ % new metric value for the current interval. Is it better?
+ if M_i<met
+ % keep this one
+ met = M_i;
+ area = A_i;
+ perimeter = P_i;
+
+ rect = [xymin;[xymax(1),xymin(2)];xymax;[xymin(1),xymax(2)];xymin];
+ rect = rect*rot';
+ rectx = rect(:,1);
+ recty = rect(:,2);
+ end
+end
+% get the final rect
+
+% all done
+
+end % mainline end
+
+
59 utils/rasterive_mesh/rasterize.m
View
@@ -0,0 +1,59 @@
+
+files=dir('*.mat');
+
+for f=1:max(size(files));
+
+fname=files(f).name;
+
+load(fname)
+
+fprintf('%s\n',fname)
+
+X=mxDomain(:,1);
+Y=mxDomain(:,2);
+[bbx,bby,~,~]=minboundrect(X,Y);
+
+% hold on
+% plot(bbx,bby,'color','r');
+
+
+
+%bottom left most point
+ [val,i]=min(bbx);
+ xll=val;
+ yll=bby(i);
+% plot(xll,yll,'o');
+%bottom right most pt
+ [val,i]=max(bbx);
+ xlr = val;
+ ylr = bby(i);
+% plot(xlr,ylr,'o');
+%top right most pt
+ [val,i]=max(bby);
+ yur = val;
+ xur = xlr;
+% plot(xur,yur,'o');
+%top left most pt
+ [val,i]=max(bby);
+ yul = val;
+ xul = xll;
+% plot(xul,yul,'o');
+
+cell_size = 2;
+
+num_rows= ceil( (yul-yll)/cell_size)
+num_cols= ceil( abs(xll-xlr)/cell_size)
+
+[mx,my]=meshgrid(xll:cell_size:xlr,yll:cell_size:yul);
+Z=gridtrimesh(tri,mxDomain,mx,my,shadows);
+
+Z=bwmorph(Z,'majority','inf');
+% imwrite(flipud(Z),strcat('model_',fname(1:end-4)),'tif');
+
+% header=[num_cols+1;num_rows+1;xll;yll;cell_size;-9999];
+% SaveAsciiRaster(flipud(Z),header,strcat('model_',fname(1:end-4)));
+end
+% imshow(Z)
+%
+
+
Please sign in to comment.
Something went wrong with that request. Please try again.