spectralpython.github.io/algorithms.html at master · spectralpython/spectralpython.github.io

History

720 lines (660 loc) · 64.7 KB

Raw

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228

229

230

231

232

233

234

235

236

237

238

239

240

241

242

243

244

245

246

247

248

249

250

251

252

253

254

255

256

257

258

259

260

261

262

263

264

265

266

267

268

269

270

271

272

273

274

275

276

277

278

279

280

281

282

283

284

285

286

287

288

289

290

291

292

293

294

295

296

297

298

299

300

301

302

303

304

305

306

307

308

309

310

311

312

313

314

315

316

317

318

319

320

321

322

323

324

325

326

327

328

329

330

331

332

333

334

335

336

337

338

339

340

341

342

343

344

345

346

347

348

349

350

351

352

353

354

355

356

357

358

359

360

361

362

363

364

365

366

367

368

369

370

371

372

373

374

375

376

377

378

379

380

381

382

383

384

385

386

387

388

389

390

391

392

393

394

395

396

397

398

399

400

401

402

403

404

405

406

407

408

409

410

411

412

413

414

415

416

417

418

419

420

421

422

423

424

425

426

427

428

429

430

431

432

433

434

435

436

437

438

439

440

441

442

443

444

445

446

447

448

449

450

451

452

453

454

455

456

457

458

459

460

461

462

463

464

465

466

467

468

469

470

471

472

473

474

475

476

477

478

479

480

481

482

483

484

485

486

487

488

489

490

491

492

493

494

495

496

497

498

499

500

501

502

503

504

505

506

507

508

509

510

511

512

513

514

515

516

517

518

519

520

521

522

523

524

525

526

527

528

529

530

531

532

533

534

535

536

537

538

539

540

541

542

543

544

545

546

547

548

549

550

551

552

553

554

555

556

557

558

559

560

561

562

563

564

565

566

567

568

569

570

571

572

573

574

575

576

577

578

579

580

581

582

583

584

585

586

587

588

589

590

591

592

593

594

595

596

597

598

599

600

601

602

603

604

605

606

607

608

609

610

611

612

613

614

615

616

617

618

619

620

621

622

623

624

625

626

627

628

629

630

631

632

633

634

635

636

637

638

639

640

641

642

643

644

645

646

647

648

649

650

651

652

653

654

655

656

657

658

659

660

661

662

663

664

665

666

667

668

669

670

671

672

673

674

675

676

677

678

679

680

681

682

683

684

685

686

687

688

689

690

691

692

693

694

695

696

697

698

699

700

701

702

703

704

705

706

707

708

709

710

711

712

713

714

715

716

717

718

719

720

<!DOCTYPE html>

<head>

<title>Spectral Algorithms — Spectral Python 0.21 documentation</title>

</head><body>

<h3>Navigation</h3>

<ul>

<a href="py-modindex.html" title="Python Module Index"

>modules</a></li>

<a href="libraries.html" title="Spectral Libraries"

accesskey="N">next</a> |</li>

<a href="graphics.html" title="Displaying Data"

accesskey="P">previous</a> |</li>

<li class="nav-item nav-item-0"><a href="index.html">Spectral Python 0.21 documentation</a> »</li>

<li class="nav-item nav-item-1"><a href="user_guide.html" accesskey="U">Spectral Python (SPy) User Guide</a> »</li>

</ul>

</div>

</div>

<h1>Spectral Algorithms<a class="headerlink" href="#spectral-algorithms" title="Permalink to this headline">¶</a></h1>

SPy implements various algorithms for dimensionality reduction and supervised &

unsupervised classification. Some of these algorithms are computationally

burdensome and require iterative access to image data. These algorithms will

almost always execute significantly faster if the image data is loaded into

memory. Therefore, if you use these algorithms and have sufficient computer

memory, you should load your image into memory.

<div class="highlight-ipython notranslate"><div class="highlight"><pre>In [1]: from spectral import *

In [2]: img = open_image('92AV3C.lan').load()

</pre></div>

</div>

<h2>Unsupervised Classification<a class="headerlink" href="#unsupervised-classification" title="Permalink to this headline">¶</a></h2>

Unsupervised classification algorithms divide image pixels into groups based on

spectral similarity of the pixels without using any prior knowledge of the

spectral classes.

<h3>k-means Clustering<a class="headerlink" href="#k-means-clustering" title="Permalink to this headline">¶</a></h3>

The k-means algorithm takes an iterative approach to generating clusters.

The parameter k specifies the desired number of clusters to generate. The algorithm

begins with an initial set of cluster centers (e.g., results from <code class="xref py py-func docutils literal notranslate">cluster</code>).

Each pixel in the image is then assigned to the nearest cluster center (using

distance in N-space as the distance metric) and each cluster center is then recomputed

as the centroid of all pixels assigned to the cluster. This process repeats until

a desired stopping criterion is reached (e.g., max number of iterations). To run

the k-means algorithm on the image and create 20 clusters, using a maximum of

50 iterations, call <a class="reference internal" href="class_func_ref.html#spectral.kmeans" title="spectral.kmeans"><code class="xref py py-func docutils literal notranslate">kmeans</code></a> as follows.

<div class="highlight-ipython notranslate"><div class="highlight"><pre>In [3]: (m, c) = kmeans(img, 20, 30)

Initializing clusters along diagonal of N-dimensional bounding box.

Starting iterations.

Iteration 1...done

21024 pixels reassigned.

Iteration 2...done

11214 pixels reassigned.

Iteration 3...done

4726 pixels reassigned.

---// snip //---

Iteration 28...done

248 pixels reassigned.

Iteration 29...done

215 pixels reassigned.

Iteration 30...done

241 pixels reassigned.

^CIteration 31... 15.9%KeyboardInterrupt: Returning clusters from previous iteration

</pre></div>

</div>

Note that I interrupted the algorithm with a keyboard interrupt (CTRL-C) after

the 30th iteration since there were only about a hundred pixels migrating between

clusters at that point. <a class="reference internal" href="class_func_ref.html#spectral.kmeans" title="spectral.kmeans"><code class="xref py py-func docutils literal notranslate">kmeans</code></a> catches the <code class="xref py py-exc docutils literal notranslate">KeyboardInterrupt</code>

exception and returns the clusters generated at the end of the previous

iteration. If you are running the algorithm interactively, this feature allows

you to set the max number of iterations to an arbitrarily high number and then

stop the algorithm when the clusters have converged to an acceptable level. If you

happen to set the max number of iterations too small (many pixels are still migrating

at the end of the final iteration), you can simply call <a class="reference internal" href="class_func_ref.html#spectral.kmeans" title="spectral.kmeans"><code class="xref py py-func docutils literal notranslate">kmeans</code></a> again

to resume processing by passing the cluster centers generated by the previous call

as the optional <code class="xref py py-class docutils literal notranslate">start_clusters</code> argument to the function.

k-means clustering results<a class="headerlink" href="#id3" title="Permalink to this image">¶</a>

</div>

Notice that <a class="reference internal" href="class_func_ref.html#spectral.kmeans" title="spectral.kmeans"><code class="xref py py-func docutils literal notranslate">kmeans</code></a> appears to have captured much more of the

spectral variation in the image than the single-pass <code class="xref py py-func docutils literal notranslate">cluster</code>

function. Let’s also take a look at the the cluster centers produced by the

algorithm:

<div class="highlight-ipython notranslate"><div class="highlight"><pre>In [4]: import matplotlib.pyplot as plt

In [5]: plt.figure()

In [6]: for i in range(c.shape[0]):

...: plt.plot(c[i])

...:

In [7]: plt.grid()

</pre></div>

</div>

k-means cluster centers<a class="headerlink" href="#id4" title="Permalink to this image">¶</a>

</div>

<h2>Supervised Classification<a class="headerlink" href="#supervised-classification" title="Permalink to this headline">¶</a></h2>

<h3>Training Data<a class="headerlink" href="#training-data" title="Permalink to this headline">¶</a></h3>

Performing supervised classification requires training a classifier with training

data that associates samples with particular training classes. To assign class

labels to pixels in an image having M rows and N columns, you must provide an

MxN integer-valued ground truth array whose elements are indices for the corresponding

training classes. A value of 0 in the ground truth array indicate an unlabeled pixel

(the pixel is not associated with a training class).

The following commands will load and display the ground truth map for our sample

image:

<div class="highlight-ipython notranslate"><div class="highlight"><pre>In [8]: gt = open_image('92AV3GT.GIS').read_band(0)

In [9]: v = imshow(classes=gt)

</pre></div>

</div>

We can now create a <code class="xref py py-class docutils literal notranslate">TrainingClassSet</code> object by calling

<a class="reference internal" href="class_func_ref.html#spectral.create_training_classes" title="spectral.create_training_classes"><code class="xref py py-func docutils literal notranslate">create_training_classes</code></a>:

<div class="highlight-ipython notranslate"><div class="highlight"><pre>In [10]: classes = create_training_classes(img, gt)

</pre></div>

</div>

With the training classes defined, we can then create a supervised classifier,

train it with the training classes, and then classify an image.

</div>

<h3>Gaussian Maximum Likelihood Classification<a class="headerlink" href="#gaussian-maximum-likelihood-classification" title="Permalink to this headline">¶</a></h3>

In this case, we’ll perform Gaussian Maximum Likelihood Classification (GMLC),

so let’s create the appropriate classifier.

<div class="highlight-ipython notranslate"><div class="highlight"><pre>In [11]: gmlc = GaussianClassifier(classes)

</pre></div>

</div>

When we created the classifier, it was automatically trained on the training sets

we provided. Notice that the classifier ignored five of the training classes.

GMLC requires computing the inverse of the covariance matrix for each training

class. Since our sample image contains 220 spectral bands, classes with fewer

than 220 samples will have singular covariance matrices, for which we can’t compute

the inverse.

Once the classifier is trained, we can use it to classify an image having

the same spectral bands as the training set. Let’s classify our training image

and display the resulting classification map.

<div class="highlight-ipython notranslate"><div class="highlight"><pre>In [12]: clmap = gmlc.classify_image(img)

In [13]: v = imshow(classes=clmap)

</pre></div>

</div>

The classification map above shows classification results for the entire image. To view

results for only the ground truth pixels we must mask out all the pixels not

associated with a training class.

<div class="highlight-ipython notranslate"><div class="highlight"><pre>In [14]: gtresults = clmap * (gt != 0)

In [15]: v = imshow(classes=gtresults)

</pre></div>

</div>

If the classification results are good, we expect the classification map

above to look very similar to the original ground truth map. To view only the

errors, we must mask out all elements in <code class="xref py py-obj docutils literal notranslate">gtResults</code> that do not match the

ground truth image.

<div class="highlight-ipython notranslate"><div class="highlight"><pre>In [16]: gterrors = gtresults * (gtresults != gt)

In [17]: v = imshow(classes=gterrors)

</pre></div>

</div>

The five contiguous regions in the error image above correspond to the ground

truth classes that the <code class="xref py py-class docutils literal notranslate">GaussianClassifier</code> ignored because they

had too few samples.

The following table lists the supervised classifiers SPy currently provides.

</colgroup>

<thead>

<tr class="row-odd"><th class="head">Classifier</th>

<th class="head">Description</th>

</tr>

</thead>

<tbody>

<tr class="row-even"><td><a class="reference internal" href="class_func_ref.html#spectral.algorithms.classifiers.GaussianClassifier" title="spectral.algorithms.classifiers.GaussianClassifier"><code class="xref py py-class docutils literal notranslate">GaussianClassifier</code></a></td>

<td>Gaussian Maximum Likelihood</td>

</tr>

<tr class="row-odd"><td><a class="reference internal" href="class_func_ref.html#spectral.algorithms.classifiers.MahalanobisDistanceClassifier" title="spectral.algorithms.classifiers.MahalanobisDistanceClassifier"><code class="xref py py-class docutils literal notranslate">MahalanobisDistanceClassifier</code></a></td>

<td>Mahalanobis Distance</td>

</tr>

<tr class="row-even"><td><code class="xref py py-class docutils literal notranslate">PerceptronClassifier</code></td>

<td>Multi-Layer Perceptron</td>

</tr>

</tbody>

</table>

</div>

<h2>Dimensionality Reduction<a class="headerlink" href="#dimensionality-reduction" title="Permalink to this headline">¶</a></h2>

Processing hyperspectral images with hundreds of bands can be computationally

burdensome and classification accuracy may suffer due to the so-called “curse

of dimensionality”. To mitigate these problems, it is often desirable to reduce

the dimensionality of the data.

<h3>Principal Components<a class="headerlink" href="#principal-components" title="Permalink to this headline">¶</a></h3>

Many of the bands within hyperspectral images are often strongly correlated.

The principal components transformation represents a linear transformation of the

original image bands to a set of new, uncorrelated features. These new features

correspond to the eigenvectors of the image covariance matrix, where the associated

eigenvalue represents the variance in the direction of the eigenvector. A very

large percentage of the image variance can be captured in a relatively small

number of principal components (compared to the original number of bands) .

The SPy function <a class="reference internal" href="class_func_ref.html#spectral.principal_components" title="spectral.principal_components"><code class="xref py py-func docutils literal notranslate">principal_components</code></a> computes the principal

components of the image data and returns the mean, covariance, eigenvalues, and

eigenvectors in a <code class="xref py py-class docutils literal notranslate">PrincipalComponents</code>.

This object also contains a transform to rotate data in to the space of the

principal compenents, as well as a method to reduce the number of eigenvectors.

<div class="highlight-ipython notranslate"><div class="highlight"><pre>In [18]: pc = principal_components(img)

In [19]: v = imshow(pc.cov)

</pre></div>

</div>

In the covariance matrix display, whiter values indicate strong positive covariance,

darker values indicate strong negative covariance, and grey values indicate

covariance near zero.

To reduce dimensionality using principal components, we can sort the eigenvalues

in descending order and then retain enough eigenvalues (an corresponding eigenvectors)

to capture a desired fraction of the total image variance. We then reduce the

dimensionality of the image pixels by projecting them onto the remaining eigenvectors.

We will choose to retain a minimum of 99.9% of the total image variance.

<div class="highlight-ipython notranslate"><div class="highlight"><pre>In [20]: pc_0999 = pc.reduce(fraction=0.999)

In [21]: # How many eigenvalues are left?

In [22]: len(pc_0999.eigenvalues)

Out[22]: 32

In [23]: img_pc = pc_0999.transform(img)

In [24]: v = imshow(img_pc[:,:,:3], stretch_all=True)

</pre></div>

</div>

Now we’ll use a Gaussian maximum likelihood classifier (GMLC) for the reduced

principal components to train and classify against the training data.

<div class="highlight-ipython notranslate"><div class="highlight"><pre>In [25]: classes = create_training_classes(img_pc, gt)

In [26]: gmlc = GaussianClassifier(classes)

In [27]: clmap = gmlc.classify_image(img_pc)

In [28]: clmap_training = clmap * (gt != 0)

In [29]: v = imshow(classes=clmap_training)

</pre></div>

</div>

And the associated errors:

<div class="highlight-ipython notranslate"><div class="highlight"><pre>In [30]: training_errors = clmap_training * (clmap_training != gt)

In [31]: v = imshow(classes=training_errors)

</pre></div>

</div>

</div>

<h3>Fisher Linear Discriminant<a class="headerlink" href="#fisher-linear-discriminant" title="Permalink to this headline">¶</a></h3>

The Fisher Linear Discriminant (a.k.a., canonical discriminant) attempts to

find a set of transformed axes that maximize the ratio of the average distance

between classes to the average distance between samples within each class. <a class="reference internal" href="#richards1999" id="id1">[Richards1999]</a>

This is written as

</div>where <img class="math" src="_images/math/cef212a0d6d030d89e5a153f271db0d6bb0a3c6f.png" alt="C_b"/> is the covariance of the class centers and <img class="math" src="_images/math/6afb20d86b0ace64a91f3099474b5ec02aacefc1.png" alt="C_w"/> is the

weighted average of the covariances of each class. If <img class="math" src="_images/math/6afb20d86b0ace64a91f3099474b5ec02aacefc1.png" alt="C_w"/> is invertible,

this becomes

</div>This eigenvalue problem is solved by the <a class="reference internal" href="class_func_ref.html#spectral.linear_discriminant" title="spectral.linear_discriminant"><code class="xref py py-func docutils literal notranslate">linear_discriminant</code></a>

function, yielding <cite>C-1</cite> eigenvalues, where <cite>C</cite> is the number of classes.

<div class="highlight-ipython notranslate"><div class="highlight"><pre>In [32]: classes = create_training_classes(img, gt)

In [33]: fld = linear_discriminant(classes)

In [34]: len(fld.eigenvectors)

Out[34]: 220

</pre></div>

</div>

Let’s view the image projected onto the top 3 components of the transform:

<div class="highlight-ipython notranslate"><div class="highlight"><pre>In [35]: img_fld = fld.transform(img)

In [36]: v = imshow(img_fld[:, :, :3])

</pre></div>

</div>

Next, we’ll classify the data using this discriminant.

<div class="highlight-ipython notranslate"><div class="highlight"><pre>In [37]: classes.transform(fld.transform)

In [38]: gmlc = GaussianClassifier(classes)

In [39]: clmap = gmlc.classify_image(img_fld)

In [40]: clmap_training = clmap * (gt != 0)

In [41]: v = imshow(classes=clmap_training)

</pre></div>

</div>

<div class="highlight-ipython notranslate"><div class="highlight"><pre>In [42]: fld_errors = clmap_training * (clmap_training != gt)

In [43]: v = imshow(classes=fld_errors)

</pre></div>

</div>

See also

<a class="reference internal" href="class_func_ref.html#spectral.algorithms.algorithms.orthogonalize" title="spectral.algorithms.algorithms.orthogonalize"><code class="xref py py-func docutils literal notranslate">orthogonalize</code></a>:

<div>Gram-Schmidt Orthogonalization

</div></blockquote>

</div>

<h2>Target Detectors<a class="headerlink" href="#target-detectors" title="Permalink to this headline">¶</a></h2>

<h3>RX Anomaly Detector<a class="headerlink" href="#rx-anomaly-detector" title="Permalink to this headline">¶</a></h3>

The RX anomaly detector uses the squared Mahalanobis distance as a measure of

how anomalous a pixel is with respect to an assumed background. The SPy

<a class="reference internal" href="class_func_ref.html#spectral.algorithms.detectors.rx" title="spectral.algorithms.detectors.rx"><code class="xref py py-func docutils literal notranslate">rx</code></a> function computes RX scores for an

array of image pixels. The squared Mahalanobis distance is given by

</div>where <img class="math" src="_images/math/888f7c323ac0341871e867220ae2d76467d74d6e.png" alt="x"/> is the pixel spectrum, <img class="math" src="_images/math/a3565136d9107b3df64b08768fe13a9bd5a8d79f.png" alt="\mu_b"/> is the background

mean, and <img class="math" src="_images/math/51c142e155733b3efe0c5fe4628975cda12a22df.png" alt="\Sigma_b"/> is the background covariance <a class="reference internal" href="#reed-yu-1990" id="id2">[Reed_Yu_1990]</a>.

If no background statistics are passed to the <code class="xref py py-func docutils literal notranslate">rx</code> function,

background statistics will be estimated from the array of pixels for which the

RX scores are to be calculated.

<div class="highlight-ipython notranslate"><div class="highlight"><pre>In [44]: rxvals = rx(img)

</pre></div>

</div>

To declare pixels as anomalous, we need to specify a threshold RX score. For

example, we could choose all image pixels whose RX score has a probability of

less than 0.001 with respect to the background:

<div class="highlight-ipython notranslate"><div class="highlight"><pre>In [45]: from scipy.stats import chi2

In [46]: nbands = img.shape[-1]

In [47]: P = chi2.ppf(0.999, nbands)

In [48]: v = imshow(1 * (rxvals > P))

</pre></div>

</div>

Rather than specifying a threshold for anomalous pixels, one can also simply

view an image of raw RX scores, where brighter pixels are considered “more anomalous”:

<div class="highlight-ipython notranslate"><div class="highlight"><pre>In [49]: v = imshow(rxvals)

</pre></div>

</div>

For the sample image, only a few pixels are visible in the image of RX scores

because a linear color scale is used and there is a very small number of pixels with RX

scores much higher than other pixels. This is apparent from viewing the histogram

of the RX scores.

<div class="highlight-ipython notranslate"><div class="highlight"><pre>In [50]: import matplotlib.pyplot as plt

In [51]: f = plt.figure()

In [52]: h = plt.hist(rxvals.ravel(), 200, log=True)

In [53]: h = plt.grid()

</pre></div>

</div>

The outliers are not obvious in the histogram, so let’s print their values:

<div class="highlight-ipython notranslate"><div class="highlight"><pre>In [54]: print(np.sort(rxvals.ravel())[-10:])

[ 665.17211462 675.85801536 683.58190673 731.4872873 739.95211335

906.98669373 956.49972325 1703.20957949 2336.11246149 9018.65517253]

</pre></div>

</div>

To view greater detail in the RX image, we can adjust the lower and upper limits

of the image display. Since we are primarily interested in the most anomalous

pixels, we will set the black level to the 99th percentile of the RX values’

cumulative histogram and set the white point to the 99.99th percentile:

<div class="highlight-ipython notranslate"><div class="highlight"><pre>In [55]: v = imshow(rxvals, stretch=(0.99, 0.9999))

</pre></div>

</div>

We can see the new RGB data limits by inspecting the returned

<a class="reference internal" href="class_func_ref.html#spectral.graphics.spypylab.ImageView" title="spectral.graphics.spypylab.ImageView"><code class="xref py py-class docutils literal notranslate">ImageView</code></a> object:

<div class="highlight-ipython notranslate"><div class="highlight"><pre>In [56]: print(v)

ImageView object:

Display bands : [0]

Interpolation : <default>

RGB data limits :

R: [291.6938067178437, 956.4997232540622]

G: [291.6938067178437, 956.4997232540622]

B: [291.6938067178437, 956.4997232540622]

</pre></div>

</div>

Note that we could also have set the contrast stretch to explicit RX values

(vice percentile values) by specifying the bounds keyword instead of stretch.

If an image contains regions with different background materials, then the

assumption of a single mean/covariance for background pixels can reduce

performance of the RX anomaly detector. In such situations, better results

can be obtained by dynamically computing background statistics in a neighborhood

around each pixel being evaluated.

To compute local background statistics for

each pixel, the <code class="xref py py-func docutils literal notranslate">rx</code> function accepts an optional <cite>window</cite> argument, which

specifies an inner/outer window within which to calculate background statistics

for each pixel being evaluated. The outer window is the window within which

background statistics will be calculated. The inner window is a smaller window

(within the outer window) indicating an exclusion zone within which pixels are

to be ignored. The purpose of the inner window is to prevent potential anomaly/target

pixels from “polluting” background statistics.

For example, to compute RX scores with background statistics computed from a

21 <img class="math" src="_images/math/8ef1f282527041d7d5c4a46dd2c60a02a7a7c00b.png" alt="\times"/> 21 pixel window about the pixel being evaluated, with an exclusion window of

5 <img class="math" src="_images/math/8ef1f282527041d7d5c4a46dd2c60a02a7a7c00b.png" alt="\times"/> 5 pixels, the function would be called as follows:

<div class="highlight-ipython notranslate"><div class="highlight"><pre>In [57]: rxvals = rx(img, window=(5,21))

</pre></div>

</div>

While the use of a windowed background will often improve results for images containing

multiple background materials, it does so at the expense of introducing two

issues. First, the sizes of the inner and outer windows must be specified such

that the resulting covariance has full rank. That is, if <img class="math" src="_images/math/50b6f08425d9c4974e96bcde8c07b8dfbf271cb8.png" alt="w_{in}"/> and

<img class="math" src="_images/math/08d5b0516b4375a7755ccd45c82be74cbf9a46e7.png" alt="w_{out}"/> represent the pixel widths of the inner and outer windows, respectively,

then <img class="math" src="_images/math/f245781aa993881751610a53d8b3be5c0671325d.png" alt="w_{out}^2 - w_{in}^2"/> must be at least as large as the number of

image bands. Second, recomputing the estimated background covariance for each

pixel in the image makes the computational complexity of of the RX score

computation orders of magnitue greater.

As a compromise between a fixed background and recomputation of mean & covariance

for each pixel, <code class="xref py py-func docutils literal notranslate">rx</code> can be passed a global covariance estimate in addition

to the window argument. In this case, only the background mean within the window will

be recomputed for each pixel. This significanly reduces computation time

for the windowed algorithm and removes the size limitation on the window (except

that the outer window must be larger than the inner). For example, since our sample image has ground

cover classes labeled, we can compute the average covariance over those ground

cover classes and use the result as an estimate of the global covariance.

<div class="highlight-ipython notranslate"><div class="highlight"><pre>In [58]: C = cov_avg(img, gt)

In [59]: rxvals = rx(img, window=(5,21), cov=C)

</pre></div>

</div>

<h3>Matched Filter<a class="headerlink" href="#matched-filter" title="Permalink to this headline">¶</a></h3>

The matched filter is a linear detector given by the formula

</div>where <img class="math" src="_images/math/bba4dc6f1475c69c3aee7dfaf608317264ba6d48.png" alt="\mu_t"/> is the target mean, <img class="math" src="_images/math/a3565136d9107b3df64b08768fe13a9bd5a8d79f.png" alt="\mu_b"/> is the background

mean, and <img class="math" src="_images/math/51c142e155733b3efe0c5fe4628975cda12a22df.png" alt="\Sigma_b"/> is the covariance. The matched filter response is

scaled such that the response is zero when the input is equal to the background

mean and equal to one when the pixel is equal to the target mean. Like the

<a class="reference internal" href="class_func_ref.html#spectral.algorithms.detectors.matched_filter" title="spectral.algorithms.detectors.matched_filter"><code class="xref py py-func docutils literal notranslate">matched_filter</code></a> function will estimate

background statistics from the input image if no background statistics are

specified.

Let’s select the image pixel at (row, col) = (8, 88) as our target, use a

global background statistics estimate, and plot all pixels whose matched filter

scores are greater than 0.2.

<div class="highlight-ipython notranslate"><div class="highlight"><pre>In [60]: t = img[8, 88]

In [61]: mfscores = matched_filter(img, t)

In [62]: v = imshow(1 * (mfscores > 0.2))

</pre></div>

</div>

As with the <code class="xref py py-func docutils literal notranslate">rx</code> function, <code class="xref py py-func docutils literal notranslate">matched_filter</code> can be applied using

windowed background statistics (optionally with a global covariance estimate).

See also

<a class="reference internal" href="class_func_ref.html#spectral.algorithms.detectors.ace" title="spectral.algorithms.detectors.ace"><code class="xref py py-func docutils literal notranslate">ace</code></a>:

<div>Adaptive Coherence/Cosine Estimator (ACE)

</div></blockquote>

</div>

<h2>Miscellaneous Functions<a class="headerlink" href="#miscellaneous-functions" title="Permalink to this headline">¶</a></h2>

<h3>Band Resampling<a class="headerlink" href="#band-resampling" title="Permalink to this headline">¶</a></h3>

Comparing spectra measured with a particular sensor to spectra collected by a

different sensor often requires resampling spectra to a common band discretization.

Spectral bands of a single sensor may drift enough over time such that spectra

collected by the same sensor at different dates requires resampling.

For resampling purposes, SPy treats a sensor as having Gaussian spectral response

functions for each of its spectral bands. A source sensor band will contribute

to any destination band where there is overlap between the FWHM of the response

functions of the two bands. If there is an overlap, an integral is

performed over the region of overlap assuming the source band data value is

constant over its FWHM (since we do not know the true spectral load over the

source band) and the destination band has a Gaussian response function. Any

target bands that do not have an overlapping source band will produce <code class="xref py py-const docutils literal notranslate">NaN</code>

as the resampled band value. If FWHM information is not available for a sensor’s

bands, each band’s FWHM is assumed to reach half the distance the its adjacent

bands. Resampling results are better when source bands are at a higher spectral

resolution than the destination bands.

To create BandResampler object, we can either pass it a

<a class="reference internal" href="class_func_ref.html#spectral.algorithms.resampling.BandResampler" title="spectral.algorithms.resampling.BandResampler"><code class="xref py py-class docutils literal notranslate">BandResampler</code></a> object for each sensor or

a list of band centers and, optionally, a list of FWHM values. Once the

BandResampler is created, we can call it with a source sensor spectrum and it will

return the resampled spectrum.

<div class="highlight-ipython notranslate"><div class="highlight"><pre>In [63]: import spectral.io.aviris as aviris

In [64]: img1 = aviris.open('f970619t01p02_r02_sc04.a.rfl', 'f970619t01p02_r02.a.spc')

In [65]: bands2 = aviris.read_aviris_bands('92AV3C.spc')

In [66]: resample = BandResampler(img1.bands, bands2)

In [67]: x1 = img1[96, 304]

In [68]: x2 = resample(x1)

</pre></div>

</div>

The Normalized Difference Vegetation Index (NDVI) is an indicator of the presence

of vegetation. The index is commonly defined as

</div>where <cite>NIR</cite> is the reflectance in the near infrared (NIR) part of the spectrum and

<cite>RED</cite> is the reflectance of the red band. For our sample image, if we take

band 21 (607.0 nm) for our red band and band 43 (802.5 nm) for near infrared,

we get the following NDVI image.

<div class="highlight-ipython notranslate"><div class="highlight"><pre>In [69]: vi = ndvi(img, 21, 43)

In [70]: v = imshow(vi)

</pre></div>

</div>

<a class="reference internal" href="class_func_ref.html#spectral.algorithms.algorithms.ndvi" title="spectral.algorithms.algorithms.ndvi"><code class="xref py py-func docutils literal notranslate">ndvi</code></a> is a simple convenience function.

You could just as easily calculate the vegetation index yourself like this:

<div class="highlight-ipython notranslate"><div class="highlight"><pre>In [71]: red = img.read_band(21)

In [72]: nir = img.read_band(43)

In [73]: vi = (nir - red) / (nir + red)

</pre></div>

</div>

<h3>Spectral Angles<a class="headerlink" href="#spectral-angles" title="Permalink to this headline">¶</a></h3>

A spectral angle refers to the angle between to spectra in N-space. In the

absence of covariance data, spectral angles can be used for classifying data

against a set of reference spectra by selecting the reference spectrum with

which the unknown spectrum has the smallest angle.

To classify our sample image to the ground truth classes using spectral angles,

we must compute the spectral angles for each pixel with each training class mean.

This is done by the <a class="reference internal" href="class_func_ref.html#spectral.algorithms.algorithms.spectral_angles" title="spectral.algorithms.algorithms.spectral_angles"><code class="xref py py-func docutils literal notranslate">spectral_angles</code></a> function.

Before calling the function, we must first create a <cite>CxB</cite> array of training class

mean spectra, where <cite>C</cite> is the number of training classes and <cite>B</cite> is the number

of spectral bands.

<div class="highlight-ipython notranslate"><div class="highlight"><pre>In [74]: import numpy as np

In [75]: classes = create_training_classes(img, gt, True)

In [76]: means = np.zeros((len(classes), img.shape[2]), float)

In [77]: for (i, c) in enumerate(classes):

....: means[i] = c.stats.mean

....:

</pre></div>

</div>

In the code above, the <code class="xref py py-const docutils literal notranslate">True</code> parameter to

of statistics for each class. Next, we call <a class="reference internal" href="class_func_ref.html#spectral.algorithms.algorithms.spectral_angles" title="spectral.algorithms.algorithms.spectral_angles"><code class="xref py py-func docutils literal notranslate">spectral_angles</code></a>,

which returns an <cite>MxNxC</cite> array, where <cite>M</cite> and <cite>N</cite> are the number of rows and

columns in the image and there are <cite>C</cite> spectral angle for each pixel. To select

the class with the smallest angle, we call the <code class="xref py py-mod docutils literal notranslate">numpy</code> <code class="xref py py-func docutils literal notranslate">argmin</code>

function to select the index for the smallest angle corresponding to each pixel.

The <code class="docutils literal notranslate">clmap + 1</code> is used in the display command because our class IDs start at 1 (not 0).

<div class="highlight-ipython notranslate"><div class="highlight"><pre>In [78]: angles = spectral_angles(img, means)

In [79]: clmap = np.argmin(angles, 2)

In [80]: v = imshow(classes=((clmap + 1) * (gt != 0)))

</pre></div>

</div>

See also

<a class="reference internal" href="class_func_ref.html#spectral.algorithms.algorithms.msam" title="spectral.algorithms.algorithms.msam"><code class="xref py py-func docutils literal notranslate">msam</code></a>: Modified Spectral Angle Mapper

</div>

References

<dt class="label" id="richards1999"><a class="fn-backref" href="#id1">Richards1999</a></dt>

<dd>Richards, J.A. & Jia, X. Remote Sensing Digital Image Analysis: An Introduction. (Springer: Berlin, 1999).

</dd>

<dd>Reed, I.S. and Yu, X., “Adaptive multiple-band CFAR detection of an optical pattern with unknown spectral distribution,” IEEE Trans. Acoust., Speech, Signal Processing, vol. 38, pp. 1760-1770, Oct. 1990.

</dd>

</dl>

</div>

<h3><a href="index.html">Table of Contents</a></h3>

<ul>

<li><a class="reference internal" href="#">Spectral Algorithms</a><ul>

<li><a class="reference internal" href="#unsupervised-classification">Unsupervised Classification</a><ul>

<li><a class="reference internal" href="#k-means-clustering">k-means Clustering</a></li>

</ul>

</li>

<li><a class="reference internal" href="#supervised-classification">Supervised Classification</a><ul>

<li><a class="reference internal" href="#training-data">Training Data</a></li>

<li><a class="reference internal" href="#gaussian-maximum-likelihood-classification">Gaussian Maximum Likelihood Classification</a></li>

</ul>

</li>

<li><a class="reference internal" href="#dimensionality-reduction">Dimensionality Reduction</a><ul>

<li><a class="reference internal" href="#principal-components">Principal Components</a></li>

<li><a class="reference internal" href="#fisher-linear-discriminant">Fisher Linear Discriminant</a></li>

</ul>

</li>

<li><a class="reference internal" href="#target-detectors">Target Detectors</a><ul>

<li><a class="reference internal" href="#rx-anomaly-detector">RX Anomaly Detector</a></li>

<li><a class="reference internal" href="#matched-filter">Matched Filter</a></li>

</ul>

</li>

<li><a class="reference internal" href="#miscellaneous-functions">Miscellaneous Functions</a><ul>

<li><a class="reference internal" href="#band-resampling">Band Resampling</a></li>

<li><a class="reference internal" href="#spectral-angles">Spectral Angles</a></li>

</ul>

</li>

</ul>

</li>

</ul>

<h4>Previous topic</h4>

<a href="graphics.html"

title="previous chapter">Displaying Data</a>

<h4>Next topic</h4>

<a href="libraries.html"

title="next chapter">Spectral Libraries</a>

<h3 id="searchlabel">Quick search</h3>

</form>

</div>

</div>

</div>

<h3>Navigation</h3>

<ul>

<a href="py-modindex.html" title="Python Module Index"

>modules</a></li>

<a href="libraries.html" title="Spectral Libraries"

>next</a> |</li>

<a href="graphics.html" title="Displaying Data"

>previous</a> |</li>

<li class="nav-item nav-item-0"><a href="index.html">Spectral Python 0.21 documentation</a> »</li>

<li class="nav-item nav-item-1"><a href="user_guide.html" >Spectral Python (SPy) User Guide</a> »</li>

</ul>

</div>

Created using <a href="http://sphinx-doc.org/">Sphinx</a> 2.4.4.

</div>

(function(i,s,o,g,r,a,m){i['GoogleAnalyticsObject']=r;i[r]=i[r]||function(){

(i[r].q=i[r].q||[]).push(arguments)},i[r].l=1*new Date();a=s.createElement(o),

m=s.getElementsByTagName(o)[0];a.async=1;a.src=g;m.parentNode.insertBefore(a,m)

})(window,document,'script','//www.google-analytics.com/analytics.js','ga');

ga('create', 'UA-47389515-2', 'spectralpython.net');

ga('send', 'pageview');

</script>

</body>

</html>

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

FilesExpand file tree

algorithms.html

Latest commit

History

algorithms.html

File metadata and controls