Hi
I'm CC'ing the list as it seems it got dropped somewhere in our
conversation.
fre, 17 04 2009 kl. 21:59 +0200, skrev Stefan Gustavson:
> I finally managed to create that OCT file for BWDIST.
> Compiles with "mkoctfile bwdist.cc". (The file "edtfunc.c"
> needs to be in the same directory.)
It works great :-)
Please commit these to SVN.
> The OCT interface was a lot nicer than the clunky old
> MEX API, and I had very little trouble getting the job
> done, but because I am not very skilled in C++ and
> have no previous experience with coding for Octave,
> I would feel more comfortable checking in my result if
> you or someone else could have a brief look at it first.
Yeah, the OCT interface is really nice. It actually makes C++
not-so-painful to use.
I've had a look at you code and I've changed a few things.
*) You can index a Matrix just like an array except you have to use
parentheses instead of brackets. This removes the need for the
calls to 'fortran_vec'.
*) Octave performs reference counting which means you don't really
have to care about handling the memory, which is really nice.
However, when you extract a Matrix (or any other object) from an
'octave_value' you get a copy of the object. However, if you ask
for a 'const Matrix' a copy isn't necessary. So, by using 'const'
you can save memory and speed things up.
*) Octave has the functions 'error' and 'warning', so you don't need
to use 'octave_stdout' for that.
*) When you extract a Matrix (or any other object) from an
'octave_value' you should check if an error occured before you do
any actions on the Matrix.
*) The license of you files said that they were part of Octave, which
really isn't true (they will be part of Octave-Forge). I've just
replaced 'Octave' with 'This program' in the license.
I've changed these things, and I'm attaching the code, so you can verify
that I didn't screw up anything.
> The function itself is tested and should be OK, but is
> my coding standard acceptable,
Octave is more strict wrt. coding style than Octave-Forge. I prefer if
we use also use the Octave coding style in the 'image' package, but it's
not a requirement. If you feel like changing the style, then it would be
nice, but if you don't want to spend time on it, then that's acceptable.
> and did I get the tex-info documentation header right?
Almost :-) You should use @var around variable names when you talk about
them in the help text also. I've changed this.
> Also, I am not 100% sure I handled the memory
> allocation correctly. Am I supposed to do any cleanup
> before exit, or is that magically taken care of by using
> OCTAVE_LOCAL_BUFFER, as I assumed?
You handled it correctly. OCTAVE_LOCAL_BUFFER cleans up after itself.
Objects like Matrix etc. are reference counted so, you don't really need
to think of them either. The OCT interface is really nice wrt. memory
handling :-)
> The comments in the source now name me as the
> copyright holder, but if common practice is to reassign
> copyright to the project manager, I would be happy
> to do that as well.
We don't do copyright assignment, so no need for this. You should have
copyright of your code.
Søren
// bwdist.cpp - OCT file, implements the BWDIST function
// Depends on "edtfunc.c" for the actual computations
/*
Copyright (C) 2009 Stefan Gustavson ([email protected])
This program is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 3 of the License, or (at your
option) any later version.
This program is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with Octave; see the file COPYING. If not, see
<http://www.gnu.org/licenses/>.
*/
#include <octave/oct.h>
#ifdef __cplusplus
extern "C"
{
#endif
#define DIST_EUCLIDEAN(x,y) ((x)*(x) + (y)*(y))
#define MAX(x,y) ((x)>(y) ? (x) : (y))
#define DIST_CHESSBOARD(x,y) (MAX(abs(x), abs(y)))
#define DIST_CITYBLOCK(x,y) (abs(x) + abs(y))
#define SQRT2_1 0.4142136536
#define DIST_QUASI_EUCLIDEAN(x,y) (abs(x)>abs(y) ? (abs(x) + SQRT2_1 * abs(y)) : (SQRT2_1 * abs(x) + abs(y)))
#define DIST(x,y) DIST_EUCLIDEAN(x,y)
#define FUNCNAME euclidean
#include "edtfunc.c"
#undef DIST
#undef FUNCNAME
#define DIST(x,y) DIST_CHESSBOARD(x,y)
#define FUNCNAME chessboard
#include "edtfunc.c"
#undef DIST
#undef FUNCNAME
#define DIST(x,y) DIST_CITYBLOCK(x,y)
#define FUNCNAME cityblock
#include "edtfunc.c"
#undef DIST
#undef FUNCNAME
#define DIST(x,y) DIST_QUASI_EUCLIDEAN(x,y)
#define FUNCNAME quasi_euclidean
#include "edtfunc.c"
#undef DIST
#undef FUNCNAME
#ifdef __cplusplus
} /* end extern "C" */
#endif
DEFUN_DLD (bwdist, args, nargout,
"-*- texinfo -*-\n\
@deftypefn {Function File} {...@var{d} =} bwdist(@var{bw})\n\
Computes the distance transform of the image @var{bw}.\n\
@var{bw} should be a binary 2D array, either a Boolean array or a\n\
numeric array containing only the values 0 and 1.\n\
The return value @var{D} is a double matrix of the same size as @var{bw}.\n\
Elements with value 0 are considered background pixels, elements\n\
with value 1 are considered object pixels. The return value\n\
for each background pixel is the distance (according to the chosen\n\
metric) to the closest object pixel. For each object pixel the\n\
return value is 0.\n\
\n\
@deftypefnx{Function File} {...@var{d} =} bwdist(@var{bw}, @var{method})\n\
\n\
@var{method} is a string to choose the distance metric. Currently\n\
available metrics are 'euclidean', 'chessboard', 'cityblock' and\n\
'quasi-euclidean', which may each be abbreviated\n\
to any string starting with 'e', 'ch', 'ci' and 'q', respectively.\n\
If @var{method} is not specified, 'euclidean' is the default.\n\
\n\
@deftypefnx {Function File} {...@var{d},@var{C}] =} bwdist(@var{bw}, @var{method})\n\
\n\
If a second output argument is given, the linear index for the\n\
closest object pixel is returned for each pixel. (For object\n\
pixels, the index points to the pixel itself.) The return value C\n\
is a matrix the same size as @var{bw}.\n\n\
@end deftypefn")
{
const int nargin = args.length();
octave_value_list retval;
/* Check for proper number of input and output arguments */
if ((nargin < 1) || (nargin>2)) {
error ("bwdist accepts only one or two input parameters.");
}
else if (nargout > 2) {
error ("bwdist returns at most 2 output parameters.");
}
else {
/* Everything is OK to proceed */
const Matrix bw = args(0).matrix_value();
if (error_state) {
error ("bwdist: input argument must be a matrix");
return retval;
}
if(bw.any_element_not_one_or_zero()) {
warning ("bwdist input contains values other than 1 and 0.");
}
dim_vector dims = bw.dims();
int rows = dims(0);
int cols = dims(1);
int caseMethod = 0; // Default 0 means Euclidean
if(nargin > 1) {
charMatrix method = args(1).char_matrix_value();
if(method(0) == 'e') caseMethod = 0; // Euclidean;
else if (method(0) == 'c') {
if(method(1) == 'h') caseMethod = 1; // chessboard
else if(method(1) == 'i') caseMethod = 2; // cityblock
}
else if(method(0) == 'q') caseMethod = 3; // quasi-Euclidean
else {
warning ("unknown metric, using 'euclidean'");
caseMethod = 0;
}
}
if (!error_state) {
/* Allocate two arrays for temporary output values */
OCTAVE_LOCAL_BUFFER (short, xdist, dims.numel());
OCTAVE_LOCAL_BUFFER (short, ydist, dims.numel());
/* Create final output array */
Matrix D (rows, cols);
/* Call the appropriate C subroutine and compute output */
switch(caseMethod) {
case 1:
chessboard(bw, rows, cols, xdist, ydist);
for(int i=0; i<rows*cols; i++) {
D(i) = DIST_CHESSBOARD(xdist[i], ydist[i]);
}
break;
case 2:
cityblock(bw, rows, cols, xdist, ydist);
for(int i=0; i<rows*cols; i++) {
D(i) = DIST_CITYBLOCK(xdist[i], ydist[i]);
}
break;
case 3:
quasi_euclidean(bw, rows, cols, xdist, ydist);
for(int i=0; i<rows*cols; i++) {
D(i) = DIST_QUASI_EUCLIDEAN(xdist[i], ydist[i]);
}
break;
case 0:
default:
euclidean(bw, rows, cols, xdist, ydist);
/* Remember sqrt() for the final output */
for(int i=0; i<rows*cols; i++) {
D(i) = sqrt((double)DIST_EUCLIDEAN(xdist[i], ydist[i]));
}
break;
}
retval(0) = D;
if(nargout > 1) {
/* Create a second output array */
Matrix C (rows, cols);
/* Compute optional 'index to closest object pixel' */
for(int i=0; i<rows*cols; i++) {
C (i) = i+1 - xdist[i] - ydist[i]*rows;
}
retval(1) = C;
}
}
}
return retval;
}
// edtfunc.c - Euclidean distance transform of a binary image
/*
Copyright (C) 2009 Stefan Gustavson ([email protected])
This program is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 3 of the License, or (at your
option) any later version.
This program is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with Octave; see the file COPYING. If not, see
<http://www.gnu.org/licenses/>.
*/
/*
* This is a sweep-and-update Euclidean distance transform of
* a binary image. All positive pixels are considered object
* pixels, zero or negative pixels are treated as background.
*
* By Stefan Gustavson ([email protected]).
*
* Originally written in 1994, based on paper-only descriptions
* of the SSED8 algorithm, invented by Per-Erik Danielsson
* and improved by Ingemar Ragnemalm. This is a classic algorithm
* with roots in the 1980s, still very good for the 2D case.
*
* Updated in 2004 to treat pixels at image edges correctly,
* and to improve code readability.
*
* Edited in 2009 to form the foundation for Octave BWDIST:
* added #define-configurable distance measure and function name
*
*/
#ifndef DIST
#define DIST(x,y) ((x) * (x) + (y) * (y))
#endif
#ifndef FUNCNAME
#define FUNCNAME edt
#endif
void FUNCNAME(const Matrix &img, int w, int h, short *distx, short *disty)
{
int x, y, i;
int offset_u, offset_ur, offset_r, offset_rd,
offset_d, offset_dl, offset_l, offset_lu;
double olddist2, newdist2, newdistx, newdisty;
int changed;
/* Initialize index offsets for the current image width */
offset_u = -w;
offset_ur = -w+1;
offset_r = 1;
offset_rd = w+1;
offset_d = w;
offset_dl = w-1;
offset_l = -1;
offset_lu = -w-1;
/* Initialize the distance images to be all large values */
for(i=0; i<w*h; i++)
if(img(i) <= 0.0)
{
distx[i] = 20000; // Large but still representable
disty[i] = 20000; // in a short (16-bit) integer
}
else
{
distx[i] = 0;
disty[i] = 0;
}
/* Perform the transformation */
do
{
changed = 0;
/* Scan rows, except first row */
for(y=1; y<h; y++)
{
/* move index to leftmost pixel of current row */
i = y*w;
/* scan right, propagate distances from above & left */
/* Leftmost pixel is special, has no left neighbors */
olddist2 = DIST(distx[i], disty[i]);
if(olddist2 > 0) // If not already zero distance
{
newdistx = distx[i+offset_u];
newdisty = disty[i+offset_u]+1;
newdist2 = DIST(newdistx, newdisty);
if(newdist2 < olddist2)
{
distx[i]=newdistx;
disty[i]=newdisty;
olddist2=newdist2;
changed = 1;
}
newdistx = distx[i+offset_ur]-1;
newdisty = disty[i+offset_ur]+1;
newdist2 = DIST(newdistx, newdisty);
if(newdist2 < olddist2)
{
distx[i]=newdistx;
disty[i]=newdisty;
changed = 1;
}
}
i++;
/* Middle pixels have all neighbors */
for(x=2; x<w-1; x++, i++)
{
olddist2 = DIST(distx[i], disty[i]);
if(olddist2 == 0) continue; // Already zero distance
newdistx = distx[i+offset_l]+1;
newdisty = disty[i+offset_l];
newdist2 = DIST(newdistx, newdisty);
if(newdist2 < olddist2)
{
distx[i]=newdistx;
disty[i]=newdisty;
olddist2=newdist2;
changed = 1;
}
newdistx = distx[i+offset_lu]+1;
newdisty = disty[i+offset_lu]+1;
newdist2 = DIST(newdistx, newdisty);
if(newdist2 < olddist2)
{
distx[i]=newdistx;
disty[i]=newdisty;
olddist2=newdist2;
changed = 1;
}
newdistx = distx[i+offset_u];
newdisty = disty[i+offset_u]+1;
newdist2 = DIST(newdistx, newdisty);
if(newdist2 < olddist2)
{
distx[i]=newdistx;
disty[i]=newdisty;
olddist2=newdist2;
changed = 1;
}
newdistx = distx[i+offset_ur]-1;
newdisty = disty[i+offset_ur]+1;
newdist2 = DIST(newdistx, newdisty);
if(newdist2 < olddist2)
{
distx[i]=newdistx;
disty[i]=newdisty;
changed = 1;
}
}
/* Rightmost pixel of row is special, has no right neighbors */
i++;
olddist2 = DIST(distx[i], disty[i]);
if(olddist2 > 0) // If not already zero distance
{
newdistx = distx[i+offset_l]+1;
newdisty = disty[i+offset_l];
newdist2 = DIST(newdistx, newdisty);
if(newdist2 < olddist2)
{
distx[i]=newdistx;
disty[i]=newdisty;
olddist2=newdist2;
changed = 1;
}
newdistx = distx[i+offset_lu]+1;
newdisty = disty[i+offset_lu]+1;
newdist2 = DIST(newdistx, newdisty);
if(newdist2 < olddist2)
{
distx[i]=newdistx;
disty[i]=newdisty;
olddist2=newdist2;
changed = 1;
}
newdistx = distx[i+offset_u];
newdisty = disty[i+offset_u]+1;
newdist2 = DIST(newdistx, newdisty);
if(newdist2 < olddist2)
{
distx[i]=newdistx;
disty[i]=newdisty;
olddist2=newdist2;
changed = 1;
}
}
/* Move index to second rightmost pixel of current row. */
/* Rightmost pixel is skipped, it has no right neighbor. */
i = y*w + w-2;
/* scan left, propagate distance from right */
for(x=w-2; x>=0; x--, i--)
{
olddist2 = DIST(distx[i], disty[i]);
if(olddist2 == 0) continue; // Already zero distance
newdistx = distx[i+offset_r]-1;
newdisty = disty[i+offset_r];
newdist2 = DIST(newdistx, newdisty);
if(newdist2 < olddist2)
{
distx[i]=newdistx;
disty[i]=newdisty;
changed = 1;
}
}
}
/* Scan rows in reverse order, except last row */
for(y=h-2; y>=0; y--)
{
/* move index to rightmost pixel of current row */
i = y*w + w-1;
/* Scan left, propagate distances from below & right */
/* Rightmost pixel is special, has no right neighbors */
olddist2 = DIST(distx[i], disty[i]);
if(olddist2 > 0) // If not already zero distance
{
newdistx = distx[i+offset_d];
newdisty = disty[i+offset_d]-1;
newdist2 = DIST(newdistx, newdisty);
if(newdist2 < olddist2)
{
distx[i]=newdistx;
disty[i]=newdisty;
olddist2=newdist2;
changed = 1;
}
newdistx = distx[i+offset_dl]+1;
newdisty = disty[i+offset_dl]-1;
newdist2 = DIST(newdistx, newdisty);
if(newdist2 < olddist2)
{
distx[i]=newdistx;
disty[i]=newdisty;
changed = 1;
}
}
i--;
/* Middle pixels have all neighbors */
for(x=w-2; x>0; x--, i--)
{
olddist2 = DIST(distx[i], disty[i]);
if(olddist2 == 0) continue; // Already zero distance
newdistx = distx[i+offset_r]-1;
newdisty = disty[i+offset_r];
newdist2 = DIST(newdistx, newdisty);
if(newdist2 < olddist2)
{
distx[i]=newdistx;
disty[i]=newdisty;
olddist2=newdist2;
changed = 1;
}
newdistx = distx[i+offset_rd]-1;
newdisty = disty[i+offset_rd]-1;
newdist2 = DIST(newdistx, newdisty);
if(newdist2 < olddist2)
{
distx[i]=newdistx;
disty[i]=newdisty;
olddist2=newdist2;
changed = 1;
}
newdistx = distx[i+offset_d];
newdisty = disty[i+offset_d]-1;
newdist2 = DIST(newdistx, newdisty);
if(newdist2 < olddist2)
{
distx[i]=newdistx;
disty[i]=newdisty;
olddist2=newdist2;
changed = 1;
}
newdistx = distx[i+offset_dl]+1;
newdisty = disty[i+offset_dl]-1;
newdist2 = DIST(newdistx, newdisty);
if(newdist2 < olddist2)
{
distx[i]=newdistx;
disty[i]=newdisty;
changed = 1;
}
}
/* Leftmost pixel is special, has no left neighbors */
olddist2 = DIST(distx[i], disty[i]);
if(olddist2 > 0) // If not already zero distance
{
newdistx = distx[i+offset_r]-1;
newdisty = disty[i+offset_r];
newdist2 = DIST(newdistx, newdisty);
if(newdist2 < olddist2)
{
distx[i]=newdistx;
disty[i]=newdisty;
olddist2=newdist2;
changed = 1;
}
newdistx = distx[i+offset_rd]-1;
newdisty = disty[i+offset_rd]-1;
newdist2 = DIST(newdistx, newdisty);
if(newdist2 < olddist2)
{
distx[i]=newdistx;
disty[i]=newdisty;
olddist2=newdist2;
changed = 1;
}
newdistx = distx[i+offset_d];
newdisty = disty[i+offset_d]-1;
newdist2 = DIST(newdistx, newdisty);
if(newdist2 < olddist2)
{
distx[i]=newdistx;
disty[i]=newdisty;
olddist2=newdist2;
changed = 1;
}
}
/* Move index to second leftmost pixel of current row. */
/* Leftmost pixel is skipped, it has no left neighbor. */
i = y*w + 1;
for(x=1; x<w; x++, i++)
{
/* scan right, propagate distance from left */
olddist2 = DIST(distx[i], disty[i]);
if(olddist2 == 0) continue; // Already zero distance
newdistx = distx[i+offset_l]+1;
newdisty = disty[i+offset_l];
newdist2 = DIST(newdistx, newdisty);
if(newdist2 < olddist2)
{
distx[i]=newdistx;
disty[i]=newdisty;
changed = 1;
}
}
}
}
while(changed); // Sweep until no more updates are made
/* The transformation is completed. */
}
------------------------------------------------------------------------------
Stay on top of everything new and different, both inside and
around Java (TM) technology - register by April 22, and save
$200 on the JavaOne (SM) conference, June 2-5, 2009, San Francisco.
300 plus technical and hands-on sessions. Register today.
Use priority code J9JMT32. http://p.sf.net/sfu/p
_______________________________________________
Octave-dev mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/octave-dev
